mirror of
https://github.com/paboyle/Grid.git
synced 2024-09-20 09:15:38 +01:00
304 lines
10 KiB
C++
304 lines
10 KiB
C++
/*************************************************************************************
|
|
|
|
Grid physics library, www.github.com/paboyle/Grid
|
|
|
|
Source file: ./tests/solver/Test_zMADWF_prec.cc
|
|
|
|
Copyright (C) 2015
|
|
|
|
Author: Christopher Kelly <ckelly@phys.columbia.edu>
|
|
|
|
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 */
|
|
|
|
//This test computes the zMobius approximation to the Mobius action and uses it within the MADWF context to accelerate an inversion
|
|
#include <Grid/Grid.h>
|
|
|
|
using namespace std;
|
|
using namespace Grid;
|
|
|
|
struct TestParams{
|
|
bool load_config;
|
|
std::string config_file;
|
|
|
|
double mass;
|
|
|
|
std::string outer_precon;
|
|
std::string inner_precon;
|
|
|
|
int Ls_outer;
|
|
double b_plus_c_outer;
|
|
double resid_outer;
|
|
|
|
int Ls_inner;
|
|
double b_plus_c_inner; //irrelevant for ZMobius
|
|
double resid_inner;
|
|
bool zmobius_inner;
|
|
double lambda_max; //upper bound of H_T eigenvalue range required to generate zMobius approximation
|
|
|
|
TestParams(): load_config(true), config_file("ckpoint_lat.1000"), mass(0.01),
|
|
Ls_outer(24), b_plus_c_outer(2.0), resid_outer(1e-8),
|
|
Ls_inner(12), b_plus_c_inner(1.0), resid_inner(1e-8), zmobius_inner(true), lambda_max(1.42), outer_precon("Standard"), inner_precon("Standard")
|
|
{}
|
|
|
|
void write(const std::string &file) const{
|
|
XmlWriter wr(file);
|
|
#define DOIT(A) wr.writeDefault(#A, A)
|
|
DOIT(load_config);
|
|
DOIT(config_file);
|
|
DOIT(mass);
|
|
DOIT(outer_precon);
|
|
DOIT(inner_precon);
|
|
DOIT(Ls_outer);
|
|
DOIT(b_plus_c_outer);
|
|
DOIT(resid_outer);
|
|
DOIT(Ls_inner);
|
|
DOIT(b_plus_c_inner);
|
|
DOIT(resid_inner);
|
|
DOIT(zmobius_inner);
|
|
DOIT(lambda_max);
|
|
#undef DOIT
|
|
}
|
|
void read(const std::string &file){
|
|
XmlReader rd(file);
|
|
#define DOIT(A) rd.readDefault(#A, A)
|
|
DOIT(load_config);
|
|
DOIT(config_file);
|
|
DOIT(mass);
|
|
DOIT(outer_precon);
|
|
DOIT(inner_precon);
|
|
DOIT(Ls_outer);
|
|
DOIT(b_plus_c_outer);
|
|
DOIT(resid_outer);
|
|
DOIT(Ls_inner);
|
|
DOIT(b_plus_c_inner);
|
|
DOIT(resid_inner);
|
|
DOIT(zmobius_inner);
|
|
DOIT(lambda_max);
|
|
#undef DOIT
|
|
}
|
|
};
|
|
|
|
struct RunParamsPrecStd{
|
|
typedef SchurRedBlackDiagMooeeSolve<LatticeFermionD> SchurSolverType;
|
|
|
|
template<typename Action>
|
|
using HermOpType = SchurDiagMooeeOperator<Action, LatticeFermionD>;
|
|
};
|
|
|
|
struct RunParamsPrecDiagTwo{
|
|
typedef SchurRedBlackDiagTwoSolve<LatticeFermionD> SchurSolverType;
|
|
|
|
template<typename Action>
|
|
using HermOpType = SchurDiagTwoOperator<Action, LatticeFermionD>;
|
|
};
|
|
|
|
|
|
struct CGincreaseTol : public MADWFinnerIterCallbackBase{
|
|
ConjugateGradient<LatticeFermionD> &cg_inner;
|
|
RealD outer_resid;
|
|
|
|
CGincreaseTol(ConjugateGradient<LatticeFermionD> &cg_inner,
|
|
RealD outer_resid): cg_inner(cg_inner), outer_resid(outer_resid){}
|
|
|
|
void operator()(const RealD current_resid){
|
|
std::cout << "CGincreaseTol with current residual " << current_resid << " changing inner tolerance " << cg_inner.Tolerance << " -> ";
|
|
while(cg_inner.Tolerance < current_resid) cg_inner.Tolerance *= 2;
|
|
//cg_inner.Tolerance = outer_resid/current_resid;
|
|
std::cout << cg_inner.Tolerance << std::endl;
|
|
}
|
|
};
|
|
|
|
template<typename RunParamsOuter, typename RunParamsInner>
|
|
void run(const TestParams ¶ms){
|
|
RealD bmc = 1.0; //use Shamir kernel
|
|
std::vector<ComplexD> gamma_inner;
|
|
|
|
std::cout << "Compute parameters" << std::endl;
|
|
if(params.zmobius_inner){
|
|
Approx::computeZmobiusGamma(gamma_inner, params.b_plus_c_inner, params.Ls_inner, params.b_plus_c_outer, params.Ls_outer, params.lambda_max);
|
|
}else{
|
|
Approx::zolotarev_data *zdata = Approx::higham(1.0,params.Ls_inner);
|
|
gamma_inner.resize(params.Ls_inner);
|
|
for(int s=0;s<params.Ls_inner;s++) gamma_inner[s] = zdata->gamma[s];
|
|
Approx::zolotarev_free(zdata);
|
|
}
|
|
std::cout << "gamma:\n";
|
|
for(int s=0;s<params.Ls_inner;s++) std::cout << s << " " << gamma_inner[s] << std::endl;
|
|
|
|
|
|
GridCartesian* UGrid = SpaceTimeGrid::makeFourDimGrid(
|
|
GridDefaultLatt(), GridDefaultSimd(Nd, vComplexD::Nsimd()),
|
|
GridDefaultMpi());
|
|
GridRedBlackCartesian* UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
|
|
|
|
|
GridCartesian* FGrid_outer = SpaceTimeGrid::makeFiveDimGrid(params.Ls_outer, UGrid);
|
|
GridCartesian* FGrid_inner = SpaceTimeGrid::makeFiveDimGrid(params.Ls_inner, UGrid);
|
|
|
|
GridRedBlackCartesian* FrbGrid_outer = SpaceTimeGrid::makeFiveDimRedBlackGrid(params.Ls_outer, UGrid);
|
|
GridRedBlackCartesian* FrbGrid_inner = SpaceTimeGrid::makeFiveDimRedBlackGrid(params.Ls_inner, UGrid);
|
|
|
|
|
|
std::vector<int> seeds4({1, 2, 3, 4});
|
|
std::vector<int> seeds5({5, 6, 7, 8});
|
|
|
|
GridParallelRNG RNG5_outer(FGrid_outer);
|
|
RNG5_outer.SeedFixedIntegers(seeds5);
|
|
|
|
GridParallelRNG RNG4(UGrid);
|
|
RNG4.SeedFixedIntegers(seeds4);
|
|
|
|
LatticeFermionD src4(UGrid); random(RNG4,src4);
|
|
|
|
LatticeFermionD result_outer(FGrid_outer);
|
|
result_outer = Zero();
|
|
LatticeGaugeFieldD Umu(UGrid);
|
|
|
|
if(params.load_config){
|
|
FieldMetaData header;
|
|
NerscIO::readConfiguration(Umu, header, params.config_file);
|
|
|
|
for(int i=0;i<Nd;i++){
|
|
assert(header.dimension[i] == GridDefaultLatt()[i]);
|
|
}
|
|
}else{
|
|
SU<Nc>::HotConfiguration(RNG4, Umu);
|
|
}
|
|
|
|
std::cout << GridLogMessage << "Lattice dimensions: " << GridDefaultLatt()
|
|
<< " Ls: " << params.Ls_outer << std::endl;
|
|
|
|
RealD M5 = 1.8;
|
|
|
|
RealD b_outer = (params.b_plus_c_outer + bmc)/2.;
|
|
RealD c_outer = (params.b_plus_c_outer - bmc)/2.;
|
|
|
|
RealD b_inner = (params.b_plus_c_inner + bmc)/2.;
|
|
RealD c_inner = (params.b_plus_c_inner - bmc)/2.;
|
|
|
|
MobiusFermionD D_outer(Umu, *FGrid_outer, *FrbGrid_outer, *UGrid, *UrbGrid, params.mass, M5, b_outer, c_outer);
|
|
ZMobiusFermionD D_inner(Umu, *FGrid_inner, *FrbGrid_inner, *UGrid, *UrbGrid, params.mass, M5, gamma_inner, b_inner, c_inner);
|
|
|
|
LatticeFermionD src_outer(FGrid_outer);
|
|
D_outer.ImportPhysicalFermionSource(src4,src_outer); //applies D_-
|
|
|
|
//Solve using a regular even-odd preconditioned CG for the Hermitian operator
|
|
//M y = x
|
|
//Mprec y'_o = x'_o where Mprec = Doo - Doe Dee^-1 Deo and x'_o = -Doe Dee^-1 x_e + x_o
|
|
//y_o = y'_o
|
|
|
|
//(Mprec^dag Mprec) y'_o = Mprec^dag x'_o
|
|
//y'_o = (Mprec^dag Mprec)^-1 Mprec^dag x'_o
|
|
|
|
//We can get Mprec^dag x'_o from x_o from SchurRedBlackDiagMooeeSolve::RedBlackSource
|
|
ConjugateGradient<LatticeFermionD> CG_outer(params.resid_outer, 10000);
|
|
typename RunParamsOuter::SchurSolverType SchurSolver_outer(CG_outer);
|
|
|
|
LatticeFermionD tmp_e_outer(FrbGrid_outer);
|
|
LatticeFermionD src_o_outer(FrbGrid_outer);
|
|
SchurSolver_outer.RedBlackSource(D_outer, src_outer, tmp_e_outer, src_o_outer);
|
|
|
|
LatticeFermionD result_o_outer(FrbGrid_outer);
|
|
result_o_outer = Zero();
|
|
|
|
GridStopWatch CGTimer;
|
|
|
|
typename RunParamsOuter::template HermOpType<MobiusFermionD> HermOpEO_outer(D_outer);
|
|
|
|
CGTimer.Start();
|
|
CG_outer(HermOpEO_outer, src_o_outer, result_o_outer);
|
|
CGTimer.Stop();
|
|
|
|
std::cout << GridLogMessage << "Total outer CG time : " << CGTimer.Elapsed()
|
|
<< std::endl;
|
|
|
|
CGTimer.Reset();
|
|
|
|
//Solve for y using MADWF with internal preconditioning
|
|
|
|
//typedef PauliVillarsSolverRBprec<LatticeFermionD, typename RunParamsOuter::SchurSolverType> PVtype;
|
|
//PVtype PV_outer(SchurSolver_outer);
|
|
|
|
typedef PauliVillarsSolverFourierAccel<LatticeFermionD, LatticeGaugeFieldD> PVtype;
|
|
PVtype PV_outer(Umu, CG_outer);
|
|
|
|
ConjugateGradient<LatticeFermionD> CG_inner(params.resid_inner, 10000, 0);
|
|
|
|
CGincreaseTol update(CG_inner, params.resid_outer);
|
|
|
|
typename RunParamsInner::SchurSolverType SchurSolver_inner(CG_inner);
|
|
|
|
ZeroGuesser<LatticeFermion> Guess;
|
|
MADWF<MobiusFermionD, ZMobiusFermionD, PVtype, typename RunParamsInner::SchurSolverType, ZeroGuesser<LatticeFermion> > madwf(D_outer, D_inner, PV_outer, SchurSolver_inner, Guess, params.resid_outer, 100, &update);
|
|
|
|
LatticeFermionD result_MADWF(FGrid_outer);
|
|
result_MADWF = Zero();
|
|
|
|
CGTimer.Start();
|
|
madwf(src4, result_MADWF);
|
|
CGTimer.Stop();
|
|
|
|
LatticeFermionD result_o_MADWF(FrbGrid_outer);
|
|
pickCheckerboard(Odd, result_o_MADWF, result_MADWF);
|
|
|
|
std::cout << GridLogMessage << "Total MADWF time : " << CGTimer.Elapsed()
|
|
<< std::endl;
|
|
|
|
LatticeFermionD diff = result_o_MADWF - result_o_outer;
|
|
std::cout <<GridLogMessage<< "Odd-parity MADWF result norm " << norm2(result_o_MADWF)
|
|
<< " Regular result norm " << norm2(result_o_outer)
|
|
<< " Norm of diff " << norm2(diff)<<std::endl;
|
|
|
|
|
|
//std::cout << GridLogMessage << "######## Dhop calls summary" << std::endl;
|
|
//D_outer.Report();
|
|
}
|
|
|
|
|
|
|
|
int main(int argc, char** argv) {
|
|
std::cout << "Init" << std::endl;
|
|
Grid_init(&argc, &argv);
|
|
|
|
TestParams params;
|
|
|
|
if( GridCmdOptionExists(argv,argv+argc,"--params") ){
|
|
std::string pfile = GridCmdOptionPayload(argv,argv+argc,"--params");
|
|
if(pfile == "TEMPLATE"){
|
|
params.write("params.templ");
|
|
return 0;
|
|
}else{
|
|
params.read(pfile);
|
|
}
|
|
}
|
|
|
|
if(params.outer_precon == "Standard" && params.inner_precon == "Standard" ){
|
|
run<RunParamsPrecStd, RunParamsPrecStd>(params);
|
|
}else if(params.outer_precon == "DiagTwo" && params.inner_precon == "Standard"){
|
|
run<RunParamsPrecDiagTwo, RunParamsPrecStd>(params);
|
|
}else if(params.outer_precon == "Standard" && params.inner_precon == "DiagTwo"){
|
|
run<RunParamsPrecStd, RunParamsPrecDiagTwo>(params);
|
|
}else if(params.outer_precon == "DiagTwo" && params.inner_precon == "DiagTwo"){
|
|
run<RunParamsPrecDiagTwo, RunParamsPrecDiagTwo>(params);
|
|
}else assert(0);
|
|
|
|
Grid_finalize();
|
|
}
|