From 3c82d16ed80b036b0ffbba785f59b213b9961eaf Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Thu, 3 Sep 2020 22:11:17 -0400 Subject: [PATCH] 4D multigrid --- tests/solver/Test_hw_multigrid.cc | 486 +++++++++++++++++++++++------- 1 file changed, 370 insertions(+), 116 deletions(-) diff --git a/tests/solver/Test_hw_multigrid.cc b/tests/solver/Test_hw_multigrid.cc index 4096b8cf..8c806805 100644 --- a/tests/solver/Test_hw_multigrid.cc +++ b/tests/solver/Test_hw_multigrid.cc @@ -27,7 +27,7 @@ Author: Peter Boyle /* END LEGAL */ #include #include -//#include +#include #include using namespace std; @@ -397,6 +397,8 @@ public: GridBase *Coarse4D; CartesianStencil Stencil; CoarsenedMatrix &Dw; + Vector NNsv; + int *NNs; GridBase * Grid(void) { return Coarse5D; }; // this is all the linalg routines need to know @@ -409,12 +411,40 @@ public: Coarse5D(&CoarseGrid5), Dw(_Dw), geom(CoarseGrid5._ndimension), - Stencil( &CoarseGrid5,geom.npoint,Even,geom.directions,geom.displacements,0) + Stencil( &CoarseGrid5,geom.npoint,Even,geom.directions,geom.displacements,0), + NNsv(Ls) { + int surface=4; + NNs =&NNsv[0]; + for(int s=0;s=(nbasis/2)) && (b< ((nbasis/2)+NNs[s]) ) ) ) }; public: - + void Project( CoarseVector &C ) + { + const int Nsimd = CComplex::Nsimd(); + autoView(Cv,C, AcceleratorWrite); + int nb2=nbasis/2; + int Ls = this->Ls; + for(int s=0;soSites(), Nsimd, { + int sF= sU*Ls+s; + auto tmp = coalescedRead(Cv[sF]); + for(int b=NN;b compressor; Stencil.HaloExchange(in,compressor); - autoView( in_v , in, AcceleratorRead); - autoView( out_v , out, AcceleratorWrite); typedef LatticeView Aview; - std::cout << "Dw"< AcceleratorViewContainer; - for(int p=0;poSites(); // Ls loop for2D int Ls=this->Ls; - std::cout << "Dw for2d"< AcceleratorViewContainer; + for(int p=0;poSites(), b, nbasis, Nsimd, { + + typedef decltype(coalescedRead(in_v[0])) calcVector; + typedef decltype(coalescedRead(in_v[0](0))) calcComplex; int sU = sF/Ls; int s = sF%Ls; + calcComplex res = Zero(); - calcVector nbr; - int ptype; - StencilEntry *SE; - for(int point=0;point_is_local) { - nbr = coalescedReadPermute(in_v[SE->_offset],ptype,SE->_permute); - } else { - nbr = coalescedRead(Stencil.CommBuf()[SE->_offset]); - } - acceleratorSynchronise(); + if(SE->_is_local) { + nbr = coalescedReadPermute(in_v[SE->_offset],ptype,SE->_permute); + } else { + nbr = coalescedRead(CBp[SE->_offset]); + } + acceleratorSynchronise(); - for(int bb=0;bb Aggregates; + + void PromoteFromSubspace(Aggregates &_Aggregates,CoarseVector &C,FineField &F) + { + auto FineGrid4 = _Aggregates.FineGrid; + FineField F4(FineGrid4); + CoarseVector C4(Coarse4D); + for(int s=0;sLs;s++){ + ExtractSlice(C4,C,s,0); + _Aggregates.PromoteFromSubspace(C4,F4); + InsertSlice(F4,F,s,0); + } + } + void ProjectToSubspace(Aggregates &_Aggregates,CoarseVector &C,FineField &F) + { + auto FineGrid4 = _Aggregates.FineGrid; + FineField F4(FineGrid4); + CoarseVector C4(Coarse4D); + for(int s=0;sLs;s++){ + ExtractSlice(F4,F,s,0); + _Aggregates.ProjectToSubspace (C4,F4); + InsertSlice(C4,C,s,0); + } + Project(C); + } + template + void Test(Aggregates &_Aggregates,GridBase *FineGrid, Ddwf &_Ddwf) + { + typedef Lattice FineField; + CoarseVector Cin(Coarse5D); + CoarseVector Cout(Coarse5D); + CoarseVector CFout(Coarse5D); + + FineField Fin(FineGrid); + FineField Fout(FineGrid); + + + std::vector seeds({1,2,3,4,5}); + GridParallelRNG RNG(Coarse5D); RNG.SeedFixedIntegers(seeds); + + gaussian(RNG,Cin); + PromoteFromSubspace(_Aggregates,Cin,Fin); + ProjectToSubspace(_Aggregates,Cin,Fin); + + std::cout << "************ "<M(Cin,Cout); + this->Project(Cout); + std::cout << " Cout "<Mdag(Cin,Cout); + this->Project(Cout); + std::cout << " Cout "< class SolverWrapper : public LinearFunction { private: LinearOperatorBase & _Matrix; @@ -508,7 +623,6 @@ public: } }; - // Must use a non-hermitian solver template class PVdagMLinearOperator : public LinearOperatorBase { @@ -544,7 +658,6 @@ public: } }; - RealD InverseApproximation(RealD x){ return 1.0/x; } @@ -572,45 +685,12 @@ public: Cheby(MdagMOp,tmp,out); } }; - -template class MirsSmoother : public LinearFunction -{ -public: - typedef LinearOperatorBase FineOperator; - Matrix & SmootherMatrix; - FineOperator & SmootherOperator; - RealD tol; - RealD shift; - int maxit; - - MirsSmoother(RealD _shift,RealD _tol,int _maxit,FineOperator &_SmootherOperator,Matrix &_SmootherMatrix) : - shift(_shift),tol(_tol),maxit(_maxit), - SmootherOperator(_SmootherOperator), - SmootherMatrix(_SmootherMatrix) - {}; - - void operator() (const Field &in, Field &out) - { - ZeroGuesser Guess; - ConjugateGradient CG(tol,maxit,false); - - Field src(in.Grid()); - - ShiftedMdagMLinearOperator,Field> MdagMOp(SmootherMatrix,shift); - SmootherOperator.AdjOp(in,src); - Guess(src,out); - CG(MdagMOp,src,out); - } -}; - -#define GridLogLevel std::cout << GridLogMessage < -class HDCRPreconditioner : public LinearFunction< Lattice > { +class MGPreconditioner : public LinearFunction< Lattice > { public: typedef Aggregation Aggregates; - typedef CoarsenedMatrix CoarseOperator; + typedef CoarseCayleyFermion CoarseOperator; typedef typename Aggregation::CoarseVector CoarseVector; typedef typename Aggregation::CoarseMatrix CoarseMatrix; typedef typename Aggregation::FineField FineField; @@ -619,25 +699,30 @@ public: Aggregates & _Aggregates; FineOperator & _FineOperator; - FineSmoother & _Smoother; + FineSmoother & _PreSmoother; + FineSmoother & _PostSmoother; + CoarseOperator & _CoarseOperator; CoarseSolver & _CoarseSolve; int level; void Level(int lv) {level = lv; }; - - HDCRPreconditioner(Aggregates &Agg, - FineOperator &Fine, - FineSmoother &Smoother, - CoarseSolver &CoarseSolve_) + MGPreconditioner(Aggregates &Agg, + FineOperator &Fine, + FineSmoother &PreSmoother, + FineSmoother &PostSmoother, + CoarseOperator &CoarseOperator_, + CoarseSolver &CoarseSolve_) : _Aggregates(Agg), _FineOperator(Fine), - _Smoother(Smoother), + _PreSmoother(PreSmoother), + _PostSmoother(PostSmoother), + _CoarseOperator(CoarseOperator_), _CoarseSolve(CoarseSolve_), level(1) { } virtual void operator()(const FineField &in, FineField & out) { - auto CoarseGrid = _Aggregates.CoarseGrid; + auto CoarseGrid = _CoarseOperator.Grid(); CoarseVector Csrc(CoarseGrid); CoarseVector Csol(CoarseGrid); FineField vec1(in.Grid()); @@ -646,45 +731,49 @@ public: double t; // Fine Smoother t=-usecond(); - _Smoother(in,out); + _PreSmoother(in,out); t+=usecond(); - GridLogLevel << "Smoother took "<< t/1000.0<< "ms" < block ({2,2,2,2}); - const int nbasis= 8; + const int nbasis= 24; auto clatt = GridDefaultLatt(); for(int d=0;d SubspaceOp(Dw); + // How to find criticall mass? + // WilsonFermionR Dw_null(Umu,*UGrid,*UrbGrid,-0.75); // 600 iters + // WilsonFermionR Dw_null(Umu,*UGrid,*UrbGrid,-0.80); // 800 iters + // WilsonFermionR Dw_null(Umu,*UGrid,*UrbGrid,-0.82); // 1023 iters + // WilsonFermionR Dw_null(Umu,*UGrid,*UrbGrid,-0.85); // 1428 iters + // WilsonFermionR Dw_null(Umu,*UGrid,*UrbGrid,-0.87); // 1900 iters + // WilsonFermionR Dw_null(Umu,*UGrid,*UrbGrid,-0.90); // 3900 iters + // WilsonFermionR Dw_null(Umu,*UGrid,*UrbGrid,-0.92); // 6200 iters + // WilsonFermionR Dw_null(Umu,*UGrid,*UrbGrid,-0.94); // 8882 iters + WilsonFermionR Dw_null(Umu,*UGrid,*UrbGrid,-0.95); // 9170 iters + // WilsonFermionR Dw_null(Umu,*UGrid,*UrbGrid,-0.96); // 8882 iters + // WilsonFermionR Dw_null(Umu,*UGrid,*UrbGrid,-0.97); // 8406 iters + // WilsonFermionR Dw_null(Umu,*UGrid,*UrbGrid,-0.99); // 6900 iters + // WilsonFermionR Dw_null(Umu,*UGrid,*UrbGrid,-1.01); // 6397 iters + // WilsonFermionR Dw_null(Umu,*UGrid,*UrbGrid,-1.00); // 5900 iters + MdagMLinearOperator MdagM_Dw(Dw_null); - Subspace Aggregates4D(Coarse4d,UGrid,0); - Subspace Aggregates5D(Coarse5d,FGrid,0); - - assert ( (nbasis & 0x1)==0); + std::cout< WilsonCG(1.0e-10,40000); + LatticeFermion w_src(UGrid); w_src=1.0; + LatticeFermion w_res(UGrid); + WilsonCG(MdagM_Dw,w_src,w_res); + exit(0); + */ std::cout< LinOpDw(Dw); std::cout< CoarseMdagM(c_Dwf); - ConjugateGradient CoarseCG(tol,MaxIt); - // BiCGSTAB CoarseBiCGSTAB(tol,MaxIt); - - c_res=Zero(); - CoarseCG(CoarseMdagM,c_src,c_res); - - std::cout< CoarseBiCGSTAB(tol,MaxIt); + ConjugateGradient CoarseCG(tol,MaxIt); + + BiCGSTAB FineBiCGSTAB(tol,MaxIt); + ConjugateGradient FineCG(tol,MaxIt); + + NonHermitianLinearOperator FineM(Ddwf); + MdagMLinearOperator FineMdagM(Ddwf); // M^\dag M + PVdagMLinearOperator FinePVdagM(Ddwf,Dpv);// M_{pv}^\dag M + SchurDiagMooeeOperator FineDiagMooee(Ddwf); // M_ee - Meo Moo^-1 Moe + SchurDiagOneOperator FineDiagOne(Ddwf); // 1 - M_ee^{-1} Meo Moo^{-1} Moe e + + MdagMLinearOperator CoarseMdagM(c_Dwf); + PVdagMLinearOperator CoarsePVdagM(c_Dwf,c_Dpv); + + std::cout< FineDiagMooee(Ddwf); // M_ee - Meo Moo^-1 Moe + f_res_e=Zero(); + FineCG(FineDiagMooee,f_src_e,f_res_e); + */ + std::cout< PM; PM(MdagM_Dw,w_src); + std::cout< cPM; cPM(CoarseMdagM,c_src); + + std::cout< CoarseZeroGuesser; + + NormalEquations CoarseCGNE (c_Dwf,CoarseCG,CoarseZeroGuesser); + ConjugateGradient CoarseMgridCG(0.001,1000); + ChebyshevSmoother FineSmoother(0.5,60.0,14,FineM,Ddwf); + + typedef MGPreconditioner > TwoLevelHDCR; + TwoLevelHDCR TwoLevelPrecon(Aggregates4D, + FineM, + FineSmoother, + FineSmoother, + c_Dwf, + CoarseCGNE); + + TwoLevelPrecon.Level(1); + PrecGeneralisedConjugateResidualNonHermitian l1PGCR(1.0e-8,100,FineM,TwoLevelPrecon,16,16); + l1PGCR.Level(1); + + f_res=Zero(); + + CoarseCG.Tolerance=0.02; + l1PGCR(f_src,f_res); + + Grid_finalize(); + exit(0); + + std::cout< CoarseMgridBiCGSTAB(0.01,1000,false); + BiCGSTAB FineMgridPreBiCGSTAB(0.0,24,false); + BiCGSTAB FineMgridPostBiCGSTAB(0.0,24,false); + ZeroGuesser FineZeroGuesser; + + SolverWrapper FineBiCGPreSmoother( FinePVdagM, FineMgridPreBiCGSTAB, FineZeroGuesser); + SolverWrapper FineBiCGPostSmoother( FinePVdagM, FineMgridPostBiCGSTAB, FineZeroGuesser); + SolverWrapper CoarsePVdagMSolver(CoarsePVdagM,CoarseMgridBiCGSTAB,CoarseZeroGuesser); + typedef MGPreconditioner > TwoLevelMG; + + TwoLevelMG _TwoLevelMG(Aggregates4D, + FinePVdagM, + FineBiCGPreSmoother, + FineBiCGPostSmoother, + c_Dwf, + CoarsePVdagMSolver); + _TwoLevelMG.Level(1); + + PrecGeneralisedConjugateResidualNonHermitian pvPGCR(1.0e-8,100,FinePVdagM,_TwoLevelMG,16,16); + pvPGCR.Level(1); + + f_res=Zero(); + pvPGCR(f_src,f_res); std::cout<