From 2111e7ab5fbb861e79174a0a400b87ffd12c21a9 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Fri, 6 Oct 2023 21:20:21 -0400 Subject: [PATCH] Run at physical mass --- tests/debug/Test_general_coarse_hdcg_phys.cc | 423 +++++++++++++++++++ 1 file changed, 423 insertions(+) create mode 100644 tests/debug/Test_general_coarse_hdcg_phys.cc diff --git a/tests/debug/Test_general_coarse_hdcg_phys.cc b/tests/debug/Test_general_coarse_hdcg_phys.cc new file mode 100644 index 00000000..07551aed --- /dev/null +++ b/tests/debug/Test_general_coarse_hdcg_phys.cc @@ -0,0 +1,423 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_general_coarse_hdcg.cc + + Copyright (C) 2023 + +Author: Peter Boyle + + 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 +#include +#include +//#include +#include + +using namespace std; +using namespace Grid; + +template +void SaveOperator(Coarsened &Operator,std::string file) +{ +#ifdef HAVE_LIME + emptyUserRecord record; + ScidacWriter WR(Operator.Grid()->IsBoss()); + assert(Operator._A.size()==Operator.geom.npoint); + WR.open(file); + for(int p=0;p +void LoadOperator(Coarsened &Operator,std::string file) +{ +#ifdef HAVE_LIME + emptyUserRecord record; + Grid::ScidacReader RD ; + RD.open(file); + assert(Operator._A.size()==Operator.geom.npoint); + for(int p=0;p +void SaveBasis(aggregation &Agg,std::string file) +{ +#ifdef HAVE_LIME + emptyUserRecord record; + ScidacWriter WR(Agg.FineGrid->IsBoss()); + WR.open(file); + for(int b=0;b +void LoadBasis(aggregation &Agg, std::string file) +{ +#ifdef HAVE_LIME + emptyUserRecord record; + ScidacReader RD ; + RD.open(file); + for(int b=0;b class TestSolver : public LinearFunction { +public: + TestSolver() {}; + void operator() (const Field &in, Field &out){ out = Zero(); } +}; + + +RealD InverseApproximation(RealD x){ + return 1.0/x; +} + +// Want Op in CoarsenOp to call MatPcDagMatPc +template +class HermOpAdaptor : public LinearOperatorBase +{ + LinearOperatorBase & wrapped; +public: + HermOpAdaptor(LinearOperatorBase &wrapme) : wrapped(wrapme) {}; + void Op (const Field &in, Field &out) { wrapped.HermOp(in,out); } + void HermOp(const Field &in, Field &out) { wrapped.HermOp(in,out); } + void AdjOp (const Field &in, Field &out){ wrapped.HermOp(in,out); } + void OpDiag (const Field &in, Field &out) { assert(0); } + void OpDir (const Field &in, Field &out,int dir,int disp) { assert(0); } + void OpDirAll (const Field &in, std::vector &out) { assert(0); }; + void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ assert(0); } +}; +template class ChebyshevSmoother : public LinearFunction +{ +public: + using LinearFunction::operator(); + typedef LinearOperatorBase FineOperator; + FineOperator & _SmootherOperator; + Chebyshev Cheby; + ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator) : + _SmootherOperator(SmootherOperator), + Cheby(_lo,_hi,_ord,InverseApproximation) + { + std::cout << GridLogMessage<<" Chebyshev smoother order "<<_ord<<" ["<<_lo<<","<<_hi<<"]"< seeds4({1,2,3,4}); + std::vector seeds5({5,6,7,8}); + std::vector cseeds({5,6,7,8}); + + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + GridParallelRNG CRNG(Coarse5d);CRNG.SeedFixedIntegers(cseeds); + + ///////////////////////// Configuration ///////////////////////////////// + LatticeGaugeField Umu(UGrid); + + FieldMetaData header; + std::string file("ckpoint_lat.4000"); + NerscIO::readConfiguration(Umu,header,file); + + //////////////////////// Fermion action ////////////////////////////////// + MobiusFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c); + + SchurDiagMooeeOperator HermOpEO(Ddwf); + + typedef HermOpAdaptor HermFineMatrix; + HermFineMatrix FineHermOp(HermOpEO); + + LatticeFermion result(FrbGrid); result=Zero(); + + LatticeFermion src(FrbGrid); random(RNG5,src); + + // Run power method on FineHermOp + PowerMethod PM; PM(HermOpEO,src); + + + //////////////////////////////////////////////////////////// + ///////////// Coarse basis and Little Dirac Operator /////// + //////////////////////////////////////////////////////////// + typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef LittleDiracOperator::CoarseVector CoarseVector; + + NextToNextToNextToNearestStencilGeometry5D geom(Coarse5d); + NearestStencilGeometry5D geom_nn(Coarse5d); + + // Warning: This routine calls PVdagM.Op, not PVdagM.HermOp + typedef Aggregation Subspace; + Subspace Aggregates(Coarse5d,FrbGrid,cb); + + //////////////////////////////////////////////////////////// + // Need to check about red-black grid coarsening + //////////////////////////////////////////////////////////// + LittleDiracOperator LittleDiracOp(geom,FrbGrid,Coarse5d); + + bool load=true; + if ( load ) { + LoadBasis(Aggregates,"Subspace.scidac"); + LoadOperator(LittleDiracOp,"LittleDiracOp.scidac"); + } else { + Aggregates.CreateSubspaceChebyshev(RNG5,HermOpEO,nbasis, + 95.0,0.1, + // 400,200,200 -- 48 iters + // 600,200,200 -- 38 iters, 162s + // 600,200,100 -- 38 iters, 169s + // 600,200,50 -- 88 iters. 370s + 800, + 200, + 100, + 0.0); + LittleDiracOp.CoarsenOperator(FineHermOp,Aggregates); + SaveBasis(Aggregates,"Subspace.scidac"); + SaveOperator(LittleDiracOp,"LittleDiracOp.scidac"); + } + + // Try projecting to one hop only + LittleDiracOperator LittleDiracOpProj(geom_nn,FrbGrid,Coarse5d); + LittleDiracOpProj.ProjectNearestNeighbour(0.01,LittleDiracOp); // smaller shift 0.02? n + + typedef HermitianLinearOperator HermMatrix; + HermMatrix CoarseOp (LittleDiracOp); + HermMatrix CoarseOpProj (LittleDiracOpProj); + + ////////////////////////////////////////// + // Build a coarse lanczos + ////////////////////////////////////////// + Chebyshev IRLCheby(0.2,40.0,71); // 1 iter + FunctionHermOp IRLOpCheby(IRLCheby,CoarseOp); + PlainHermOp IRLOp (CoarseOp); + int Nk=48; + int Nm=64; + int Nstop=Nk; + ImplicitlyRestartedLanczos IRL(IRLOpCheby,IRLOp,Nstop,Nk,Nm,1.0e-5,20); + + int Nconv; + std::vector eval(Nm); + std::vector evec(Nm,Coarse5d); + CoarseVector c_src(Coarse5d); c_src=1.0; + CoarseVector c_res(Coarse5d); + CoarseVector c_ref(Coarse5d); + + PowerMethod cPM; cPM(CoarseOp,c_src); + + IRL.calc(eval,evec,c_src,Nconv); + DeflatedGuesser DeflCoarseGuesser(evec,eval); + + ////////////////////////////////////////// + // Build a coarse space solver + ////////////////////////////////////////// + int maxit=20000; + ConjugateGradient CG(1.0e-8,maxit,false); + ConjugateGradient CGfine(1.0e-8,10000,false); + ZeroGuesser CoarseZeroGuesser; + + // HPDSolver HPDSolve(CoarseOp,CG,CoarseZeroGuesser); + HPDSolver HPDSolve(CoarseOp,CG,DeflCoarseGuesser); + c_res=Zero(); + HPDSolve(c_src,c_res); c_ref = c_res; + std::cout << GridLogMessage<<"src norm "< HPDSolveProj(CoarseOpProj,CG,DeflCoarseGuesser); + c_res=Zero(); + HPDSolveProj(c_src,c_res); + std::cout << GridLogMessage<<"src norm "< + CoarseSmoother(1.0,37.,8,CoarseOpProj); // just go to sloppy 0.1 convergence + // CoarseSmoother(0.1,37.,8,CoarseOpProj); // + // CoarseSmoother(0.5,37.,6,CoarseOpProj); // 8 iter 0.36s + // CoarseSmoother(0.5,37.,12,CoarseOpProj); // 8 iter, 0.55s + // CoarseSmoother(0.5,37.,8,CoarseOpProj);// 7-9 iter + // CoarseSmoother(1.0,37.,8,CoarseOpProj); // 0.4 - 0.5s solve to 0.04, 7-9 iter + // ChebyshevSmoother CoarseSmoother(0.5,36.,10,CoarseOpProj); // 311 + + //////////////////////////////////////////////////////// + // CG, Cheby mode spacing 200,200 + // Unprojected Coarse CG solve to 1e-8 : 190 iters, 4.9s + // Unprojected Coarse CG solve to 4e-2 : 33 iters, 0.8s + // Projected Coarse CG solve to 1e-8 : 100 iters, 0.36s + //////////////////////////////////////////////////////// + // CoarseSmoother(1.0,48.,8,CoarseOpProj); 48 evecs + //////////////////////////////////////////////////////// + // ADEF1 Coarse solve to 1e-8 : 44 iters, 2.34s 2.1x gain + // ADEF1 Coarse solve to 4e-2 : 7 iters, 0.4s + // HDCG 38 iters 162s + // + // CoarseSmoother(1.0,40.,8,CoarseOpProj); 48 evecs + // ADEF1 Coarse solve to 1e-8 : 37 iters, 2.0s 2.1x gain + // ADEF1 Coarse solve to 4e-2 : 6 iters, 0.36s + // HDCG 38 iters 169s + + TwoLevelADEF1defl + cADEF1(1.0e-8, 500, + CoarseOp, + CoarseSmoother, + evec,eval); + + c_res=Zero(); + cADEF1(c_src,c_res); + std::cout << GridLogMessage<<"src norm "< Smoother(10.0,100.0,10,FineHermOp); //499 + // ChebyshevSmoother Smoother(3.0,100.0,10,FineHermOp); //383 + // ChebyshevSmoother Smoother(1.0,100.0,10,FineHermOp); //328 + // std::vector los({0.5,1.0,3.0}); // 147/142/146 nbasis 1 + // std::vector los({1.0,2.0}); // Nbasis 24: 88,86 iterations + // std::vector los({2.0,4.0}); // Nbasis 32 == 52, iters + // std::vector los({2.0,4.0}); // Nbasis 40 == 36,36 iters + + // + // Turns approx 2700 iterations into 340 fine multiplies with Nbasis 40 + // Need to measure cost of coarse space. + // + // -- i) Reduce coarse residual -- 0.04 + // -- ii) Lanczos on coarse space -- done + // -- iii) Possible 1 hop project and/or preconditioning it - easy - PrecCG it and + // use a limited stencil. Reread BFM code to check on evecs / deflation strategy with prec + // + std::vector los({3.0}); // Nbasis 40 == 36,36 iters + + // std::vector ords({7,8,10}); // Nbasis 40 == 40,38,36 iters (320,342,396 mults) + std::vector ords({7}); // Nbasis 40 == 40 iters (320 mults) + + for(int l=0;l CGsloppy(4.0e-2,maxit,false); + HPDSolver HPDSolveSloppy(CoarseOp,CGsloppy,DeflCoarseGuesser); + + // ChebyshevSmoother Smoother(lo,92,10,FineHermOp); // 36 best case + ChebyshevSmoother Smoother(lo,92,ords[o],FineHermOp); // 311 + + ////////////////////////////////////////// + // Build a HDCG solver + ////////////////////////////////////////// + TwoLevelADEF2 + HDCG(1.0e-8, 100, + FineHermOp, + Smoother, + HPDSolveSloppy, + HPDSolve, + Aggregates); + + TwoLevelADEF2 + HDCGdefl(1.0e-8, 100, + FineHermOp, + Smoother, + cADEF1, + HPDSolve, + Aggregates); + + result=Zero(); + HDCGdefl(src,result); + + result=Zero(); + HDCG(src,result); + + + } + } + + // Standard CG + result=Zero(); + CGfine(HermOpEO, src, result); + + Grid_finalize(); + return 0; +}