diff --git a/lib/util/Init.cc b/lib/util/Init.cc index 1266d34d..3bdd52ad 100644 --- a/lib/util/Init.cc +++ b/lib/util/Init.cc @@ -157,6 +157,14 @@ void GridCmdOptionInt(std::string &str,int & val) return; } +// ypj [add] +void GridCmdOptionFloat(std::string &str,double & val) +{ + std::stringstream ss(str); + ss>>val; + return; +} + void GridParseLayout(char **argv,int argc, std::vector &latt, diff --git a/lib/util/Init.h b/lib/util/Init.h index 3da00742..e71c95c6 100644 --- a/lib/util/Init.h +++ b/lib/util/Init.h @@ -54,6 +54,9 @@ namespace Grid { std::string GridCmdVectorIntToString(const std::vector & vec); void GridCmdOptionCSL(std::string str,std::vector & vec); void GridCmdOptionIntVector(std::string &str,std::vector & vec); + // ypj [add] + void GridCmdOptionInt(std::string &str,int & val); + void GridCmdOptionFloat(std::string &str,double & val); void GridParseLayout(char **argv,int argc, diff --git a/tests/lanczos/Test_dwf_block_lanczos.cc b/tests/lanczos/Test_dwf_block_lanczos.cc index 64b72c98..71879d25 100644 --- a/tests/lanczos/Test_dwf_block_lanczos.cc +++ b/tests/lanczos/Test_dwf_block_lanczos.cc @@ -35,17 +35,125 @@ typedef typename GparityDomainWallFermionR::FermionField FermionField; RealD AllZero(RealD x){ return 0.;} +class CmdJobParams +{ + public: + std::string gaugefile; + + int Ls; + double mass; + double M5; + double mob_b; + + int Nu; + int Nk; + int Np; + int Nstop; + int MaxIter; + double resid; + + double low; + double high; + int order; + + CmdJobParams() + : gaugefile("Hot"), + Ls(8), mass(0.01), M5(1.8), mob_b(1.5), + Nu(4), Nk(200), Np(200), Nstop(100), MaxIter(10), resid(1.0e-8), + low(0.2), high(5.5), order(11) + {}; + + void Parse(char **argv, int argc); +}; + + +void CmdJobParams::Parse(char **argv,int argc) +{ + std::string arg; + std::vector vi; + + if( GridCmdOptionExists(argv,argv+argc,"--gconf") ){ + gaugefile = GridCmdOptionPayload(argv,argv+argc,"--gconf"); + } + + if( GridCmdOptionExists(argv,argv+argc,"--Ls") ){ + arg = GridCmdOptionPayload(argv,argv+argc,"--Ls"); + GridCmdOptionInt(arg,Ls); + } + + if( GridCmdOptionExists(argv,argv+argc,"--mass") ){ + arg = GridCmdOptionPayload(argv,argv+argc,"--mass"); + GridCmdOptionFloat(arg,mass); + } + + if( GridCmdOptionExists(argv,argv+argc,"--M5") ){ + arg = GridCmdOptionPayload(argv,argv+argc,"--M5"); + GridCmdOptionFloat(arg,M5); + } + + if( GridCmdOptionExists(argv,argv+argc,"--mob_b") ){ + arg = GridCmdOptionPayload(argv,argv+argc,"--mob_b"); + GridCmdOptionFloat(arg,mob_b); + } + + if( GridCmdOptionExists(argv,argv+argc,"--irbl") ){ + arg = GridCmdOptionPayload(argv,argv+argc,"--irbl"); + GridCmdOptionIntVector(arg,vi); + Nu = vi[0]; + Nk = vi[1]; + Np = vi[2]; + Nstop = vi[3]; + MaxIter = vi[4]; + } + + if( GridCmdOptionExists(argv,argv+argc,"--resid") ){ + arg = GridCmdOptionPayload(argv,argv+argc,"--resid"); + GridCmdOptionFloat(arg,resid); + } + + if( GridCmdOptionExists(argv,argv+argc,"--cheby_l") ){ + arg = GridCmdOptionPayload(argv,argv+argc,"--cheby_l"); + GridCmdOptionFloat(arg,low); + } + + if( GridCmdOptionExists(argv,argv+argc,"--cheby_u") ){ + arg = GridCmdOptionPayload(argv,argv+argc,"--cheby_u"); + GridCmdOptionFloat(arg,high); + } + + if( GridCmdOptionExists(argv,argv+argc,"--cheby_n") ){ + arg = GridCmdOptionPayload(argv,argv+argc,"--cheby_n"); + GridCmdOptionInt(arg,order); + } + + if ( CartesianCommunicator::RankWorld() == 0 ) { + clog <<" Gauge Configuration "<< gaugefile << '\n'; + clog <<" Ls "<< Ls << '\n'; + clog <<" mass "<< mass << '\n'; + clog <<" M5 "<< M5 << '\n'; + clog <<" mob_b "<< mob_b << '\n'; + clog <<" Nu "<< Nu << '\n'; + clog <<" Nk "<< Nk << '\n'; + clog <<" Np "<< Np << '\n'; + clog <<" Nstop "<< Nstop << '\n'; + clog <<" MaxIter "<< MaxIter << '\n'; + clog <<" resid "<< resid << '\n'; + clog <<" Cheby Poly "<< low << "," << high << "," << order << std::endl; + } +} + + int main (int argc, char ** argv) { Grid_init(&argc,&argv); - - const int Ls=8; - //const int Ls=16; + + CmdJobParams JP; + JP.Parse(argv,argc); 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); + GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(JP.Ls,UGrid); + GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(JP.Ls,UGrid); printf("UGrid=%p UrbGrid=%p FGrid=%p FrbGrid=%p\n",UGrid,UrbGrid,FGrid,FrbGrid); std::vector seeds4({1,2,3,4}); @@ -56,31 +164,23 @@ int main (int argc, char ** argv) // ypj [note] why seed RNG5 again? bug? In this case, run with a default seed(). //GridParallelRNG RNG5rb(FrbGrid); //RNG5rb.SeedFixedIntegers(seeds5); -#if 0 // ypj [note] generate a random configuration - LatticeGaugeField Umu(UGrid); - SU3::HotConfiguration(RNG4, Umu); - - std::vector U(4,UGrid); - for(int mu=0;mu(Umu,mu); - } -#else // read configuration from a file LatticeGaugeField Umu(UGrid); std::vector U(4,UGrid); - FieldMetaData header; - std::string file("./ckpoint_lat"); - NerscIO::readConfiguration(Umu,header,file); - + if ( JP.gaugefile.compare("Hot") == 0 ) { + SU3::HotConfiguration(RNG4, Umu); + } else { + FieldMetaData header; + NerscIO::readConfiguration(Umu,header,JP.gaugefile); + } + for(int mu=0;mu(Umu,mu); } -#endif - //RealD mass=0.01; - RealD mass=0.00107; - RealD M5=1.8; - RealD mob_b=1.5; + RealD mass = JP.mass; + RealD M5 = JP.M5; + RealD mob_b = JP.mob_b; // DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); GparityMobiusFermionD ::ImplParams params; std::vector twists({1,1,1,0}); @@ -92,27 +192,25 @@ int main (int argc, char ** argv) SchurDiagTwoOperator HermOp(Ddwf); // SchurDiagMooeeOperator HermOp(Ddwf); - const int Nstop = 100; - const int Nu = 10; - const int Nk = 200; - const int Np = 200; - const int Nm = Nk+Np; - const int MaxIt= 10; - RealD resid = 1.0e-8; + int Nu = JP.Nu; + int Nk = JP.Nk; + int Nm = Nk+JP.Np; //std::vector Coeffs { 0.,-1.}; // ypj [note] this may not be supported by some compilers std::vector Coeffs({ 0.,-1.}); Polynomial PolyX(Coeffs); //Chebyshev Cheb(0.2,5.5,11); - Chebyshev Cheb(0.03,5.5,21); -// ChebyshevLanczos Cheb(9.,1.,0.,20); + Chebyshev Cheb(JP.low,JP.high,JP.order); // Cheb.csv(std::cout); -// exit(-24); - ImplicitlyRestartedBlockLanczos IRBL(HermOp,Cheb,Nstop,Nu,Nk,Nm,resid,MaxIt); - + ImplicitlyRestartedBlockLanczos IRBL(HermOp, + Cheb, + JP.Nstop, + Nu,Nk,Nm, + JP.resid, + JP.MaxIter); - std::vector eval(Nm); + std::vector eval(Nm); std::vector src(Nu,FrbGrid); for ( int i=0; i