diff --git a/Grid/algorithms/multigrid/MrhsMultiGrid.h b/Grid/algorithms/multigrid/MrhsMultiGrid.h new file mode 100644 index 000000000..e8edf9f5b --- /dev/null +++ b/Grid/algorithms/multigrid/MrhsMultiGrid.h @@ -0,0 +1,256 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./Grid/algorithms/multigrid/MrhsMultiGrid.h + + Copyright (C) 2026 + +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. + + See the full license in the file "LICENSE" in the top level distribution + directory +*************************************************************************************/ +/* END LEGAL */ +#pragma once + +#include + +NAMESPACE_BEGIN(Grid); + +////////////////////////////////////////////////////////////////////// +// mrhs LinearFunction interface: vector-of-fields in, vector out. +// The outer level carries mrhs as std::vector; below it mrhs +// is PACKED into a single D+1 field (rhs = dim 0) and the coarse +// classes are plain LinearFunctions on that. +////////////////////////////////////////////////////////////////////// +template +class MrhsLinearFunction { +public: + virtual void operator()(std::vector &in, std::vector &out) = 0; +}; + +////////////////////////////////////////////////////////////////////// +// Single-polynomial mrhs PGCR: one step length and one set of +// orthogonalisation coefficients SHARED across the right-hand sides +// (vnorm2/vinnerProduct sum over rhs). For DWF the dense spectrum +// makes sharing near-free and the summed reductions amortise; see the +// mrhs-HDCG paper (arXiv:2409.03904). Orthogonalisation is classical +// Gram-Schmidt, batched per rhs (rankInnerProductMulti/axpyMulti) with +// ONE GlobalSumVector for all coefficients. +// +// OnStep is a programmatic per-step hook for algorithmic studies +// (e.g. driving the GCR coefficient recorder); it is not a consumer +// parameter. +////////////////////////////////////////////////////////////////////// +template +class MrhsPGCRNonHermitian { +public: + RealD Tolerance; Integer MaxIterations; int mmax,nstep,steps,level; + int ZeroGuess = 0; int FirstCycle = 0; + std::string name = "Level 1"; + LinearOperatorBase &Linop; + MrhsLinearFunction &Preconditioner; + std::function OnStep; // called with the outer step count after every step + void Level(int lv){ name = "Level " + std::to_string(lv); level=lv; } + void Name(std::string n){ name = n; } + void SetZeroGuess(int z){ ZeroGuess=z; } + MrhsPGCRNonHermitian(RealD tol,Integer maxit,LinearOperatorBase &_Linop,MrhsLinearFunction &Prec,int _mmax,int _nstep) + : Tolerance(tol),MaxIterations(maxit),Linop(_Linop),Preconditioner(Prec),mmax(_mmax),nstep(_nstep){ level=1; } + static RealD vnorm2(std::vector &x){ RealD s=0; for(auto &f:x) s+=norm2(f); return s; } + static ComplexD vinnerProduct(std::vector &x,std::vector &y){ ComplexD s(0); for(int r=0;r<(int)x.size();r++) s+=innerProduct(x[r],y[r]); return s; } + static void vaxpy(std::vector &z,ComplexD a,std::vector &x,std::vector &y){ for(int r=0;r<(int)z.size();r++) axpy(z[r],a,x[r],y[r]); } + void vOp(std::vector &in,std::vector &out){ GRID_TRACE("MrhsPGCR::vOp"); for(int r=0;r<(int)in.size();r++) Linop.Op(in[r],out[r]); } + void operator()(std::vector &src,std::vector &psi){ + RealD cp,ssq,rsq; int nrhs=src.size(); GridBase *grid=src[0].Grid(); + ssq=vnorm2(src); rsq=Tolerance*Tolerance*ssq; + std::vector r(nrhs,grid); + GridStopWatch T; T.Start(); steps=0; FirstCycle=1; + for(int k=0;k &src,std::vector &psi,RealD rsq){ + RealD cp; ComplexD a,rq; int nrhs=src.size(); GridBase *grid=src[0].Grid(); + std::vector r(nrhs,grid),Az(nrhs,grid); // Az: restart residual scratch only + std::vector< std::vector > q(mmax,std::vector(nrhs,grid)); + std::vector< std::vector > p(mmax,std::vector(nrhs,grid)); + std::vector qq(mmax); + if (ZeroGuess && FirstCycle) { for(int rr=0;rr(mmax-1))?(mmax-1):(kp); + { + GRID_TRACE("MrhsPGCR orthog"); + // Classical Gram-Schmidt: all coefficients against the UN-updated new q + // (independent, batchable), then apply. Complex coefficient: the + // operator is non-Hermitian, real() alone left q's non-orthogonal. + // Batched per rhs (one fused kernel + one reduction each), the shared + // coefficient summed over rhs on the host, ONE GlobalSumVector. + std::vector bcoef(northog,ComplexD(0.0)), part; + for(int rr=0;rr qwin(northog); + for(int back=0;back=0); qwin[back]=&q[peri_back][rr]; } + rankInnerProductMulti(part,qwin,q[peri_kp][rr]); + for(int back=0;backGlobalSumVector(&bcoef[0],northog); + for(int back=0;back qwin(northog), pwin(northog); + for(int back=0;back +class MrhsDenseCCSolve : public LinearFunction { +public: + DenseType &_Dense; + int _nrhs; + MrhsDenseCCSolve(DenseType &D, int nrhs) : _Dense(D), _nrhs(nrhs) {} + using LinearFunction::operator(); + virtual void operator()(const CoarseCoarseField &in, CoarseCoarseField &out){ + _Dense.ApplyBatch6D(in, out, _nrhs); + } +}; + +////////////////////////////////////////////////////////////////////// +// L2->L3 half V-cycle on the D+1 coarse field: coarse-coarse correction +// through the mixed blockProject, then post-smooth. +////////////////////////////////////////////////////////////////////// +template +class MrhsCoarseThreeLevelPrec : public LinearFunction { +public: + LinearOperatorBase &_CoarseOp; + LinearFunction &_CoarseSmoother; + MultiRHSBlockProject &_Projector; + LinearFunction &_CoarseCoarseSolve; + GridBase *_Coarse5d, *_CoarseCoarse5d, *_CoarseCoarseMrhs; + int _nrhs; + MrhsCoarseThreeLevelPrec(LinearOperatorBase &CoarseOp, + LinearFunction &CoarseSmoother, + MultiRHSBlockProject &Projector, + LinearFunction &CoarseCoarseSolve, + GridBase *Coarse5d, GridBase *CoarseCoarse5d, GridBase *CoarseCoarseMrhs, int nrhs) + : _CoarseOp(CoarseOp), _CoarseSmoother(CoarseSmoother), _Projector(Projector), + _CoarseCoarseSolve(CoarseCoarseSolve), + _Coarse5d(Coarse5d), _CoarseCoarse5d(CoarseCoarse5d), _CoarseCoarseMrhs(CoarseCoarseMrhs), _nrhs(nrhs) {} + using LinearFunction::operator(); + virtual void operator()(const CoarseField &in, CoarseField &out) { + CoarseField vec1(in.Grid()); + CoarseField vec2(in.Grid()); + out = in; + _CoarseOp.Op(out,vec1); sub(vec1,in,vec1); + + // restrict, through the mixed blockProject: D+1 coarse in, D+1 cc out + CoarseCoarseField CCsrc(_CoarseCoarseMrhs); + CoarseCoarseField CCsol(_CoarseCoarseMrhs); + _Projector.blockProject(vec1,CCsrc); + + _CoarseCoarseSolve(CCsrc,CCsol); + + _Projector.blockPromote(vec1,CCsol); + add(out,out,vec1); + + _CoarseOp.Op(out,vec1); sub(vec1,in,vec1); + _CoarseSmoother(vec1,vec2); + add(out,out,vec2); + } +}; + +////////////////////////////////////////////////////////////////////// +// L1->L2 mrhs V-cycle. The whole V-cycle is preconditioner: its fine +// residuals and the smoother may run with sloppy halos; the caller +// (the outer Krylov) gets the exact operator back on exit. SetSloppy +// is wired by the composer to PVdagMLinearOperator::SloppyComms (a +// no-op by default), replacing the file-scope global the example used. +////////////////////////////////////////////////////////////////////// +template +class MrhsTwoLevelMG : public MrhsLinearFunction { +public: + typedef MrhsCoarseVector CoarseVector; + LinearOperatorBase &_FineOperator; + FineSmoother &_PostSmoother; + MultiRHSBlockProject &_Projector; + LinearFunction &_CoarseSolve; + GridBase *_CoarseGrid, *_CoarseGridMrhs; + std::function SetSloppy = [](int){}; + int SloppyComms = 0; // value passed to SetSloppy on entry + MrhsTwoLevelMG(LinearOperatorBase &FineOp, FineSmoother &Post, + MultiRHSBlockProject &Projector, LinearFunction &CoarseSolve, + GridBase *CoarseGrid, GridBase *CoarseGridMrhs) + : _FineOperator(FineOp),_PostSmoother(Post),_Projector(Projector),_CoarseSolve(CoarseSolve), + _CoarseGrid(CoarseGrid),_CoarseGridMrhs(CoarseGridMrhs){} + virtual void operator()(std::vector &in, std::vector &out){ + GRID_TRACE("MGVcycle"); + SetSloppy(SloppyComms); + int nrhs=in.size(); GridBase *fgrid=in[0].Grid(); + std::vector vec1(nrhs,fgrid),vec2(nrhs,fgrid); + for(int r=0;r D+1 coarse, via the mixed blockProject + CoarseVector CsrcMrhs(_CoarseGridMrhs), CsolMrhs(_CoarseGridMrhs); + { GRID_TRACE("MGProject"); + _Projector.blockProject(vec1,CsrcMrhs); + } + CsolMrhs=Zero(); + { GRID_TRACE("MGCoarseSolve"); + _CoarseSolve(CsrcMrhs,CsolMrhs); + } + { GRID_TRACE("MGPromote"); + _Projector.blockPromote(vec1,CsolMrhs); + for(int r=0;r + + 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. + + See the full license in the file "LICENSE" in the top level distribution + directory +*************************************************************************************/ +/* END LEGAL */ +#pragma once + +#include +#include +#include + +NAMESPACE_BEGIN(Grid); + +////////////////////////////////////////////////////////////////////////////////////// +// The three-level mrhs PVdagM multigrid, as objects: +// +// MGCoarseGrids the derived coarse grids (owned here, since +// conformability is pointer identity and someone +// must hold them; declare it BEFORE anything that +// borrows from it) +// PVdagMMultiGridCoarsening ALL the coarsening information for one gauge +// configuration: the raw near-null basis, the +// transfer operators it induces, the Galerkin +// coarse operator at each level, and the dense +// bottom inverse. Everything here is a function +// of the gauge field -- the state HMC must +// rebuild or maintain as U evolves -- and +// everything is solver-independent. +// PVdagMMultiGridSolver the solve chain composed on a borrowed +// coarsening: smoothers, coarse Krylov, V-cycle, +// outer mrhs PGCR. +// +// Scope discipline: grids outlive the coarsening outlives the solver. +////////////////////////////////////////////////////////////////////////////////////// + +////////////////////////////////////////////////////////////////////// +// Subspace I/O: bare-vector scidac records. Loaded vectors are RAW -- +// deliberately NOT re-orthogonalised: CoarsenOperator block +// orthonormalises in place, and projecting a block-orthonormal vector +// onto its own block-orthonormalised aggregation gives e_k with the +// near-null content silently gone. GramGuard below catches that. +////////////////////////////////////////////////////////////////////// +template +void saveSubspace(std::vector &subspace, std::string const fname){ +#ifdef HAVE_LIME + Grid::emptyUserRecord record; + Grid::ScidacWriter SW(subspace[0].Grid()->IsBoss()); + SW.open(fname); + for (int k = 0; k < (int)subspace.size(); k++) { + SW.writeScidacFieldRecord(subspace[k], record); + } + SW.close(); +#endif +} +template +void loadSubspace(std::vector &subspace, std::string const fname){ +#ifdef HAVE_LIME + Grid::emptyUserRecord record; + Grid::ScidacReader SR; + SR.open(fname); + for (int k = 0; k < (int)subspace.size(); k++) { + SR.readScidacFieldRecord(subspace[k], record); + } + SR.close(); +#endif +} + +////////////////////////////////////////////////////////////////////// +// || - I||_F over a set of coarse vectors. Small means the raw +// near-null content survived the projection. +////////////////////////////////////////////////////////////////////// +template +RealD GramDefect(std::vector &v) +{ + RealD s2=0.0; + for(int i=0;i<(int)v.size();i++){ + for(int j=0;j<(int)v.size();j++){ + ComplexD sij=TensorRemove(innerProduct(v[i],v[j])); + ComplexD d=sij-(i==j?ComplexD(1.0):ComplexD(0.0)); + s2+=real(d)*real(d)+imag(d)*imag(d); + } + } + return std::sqrt(s2); +} + +// On a leak every image collapses to the block unit e_k, the Gram becomes +// N*I, and the defect lands at (N-1)*sqrt(nbasis) -- orders above the ~0.2 +// of a content-preserving projection. Trip well below that so a mis-set +// threshold costs a log line rather than the run. +template +void GramGuard(const std::string &name,std::vector &v,GridBase *grid) +{ + RealD defect = GramDefect(v); + RealD N = (RealD)grid->gSites(); + RealD leak = (N-1.0)*std::sqrt((RealD)v.size()); + RealD trip = std::sqrt(N); + std::cout << GridLogMessage << "GUARD: ||<"< - I||_F = " << defect + << " (e_k leak would be " << leak << ", trip at " << trip << ")" << std::endl; + GRID_ASSERT( defect < trip ); +} + +////////////////////////////////////////////////////////////////////// +// The derived grids of the three-level chain, owned in one place. +// The coarse space is UNVECTORISED (sComplex scalar, simd {1,..,1}); +// the 5D/6D grids are built directly so the SIMD layout is ours. +// rhs/batch is dim 0 of the 6D grids, undistributed -- no divisibility +// constraint on nrhs. +////////////////////////////////////////////////////////////////////// +class MGCoarseGrids { +public: + GridCartesian *FGrid; // borrowed + int Ls; + int batch; + Coordinate clatt; // 4d coarse lattice + Coordinate cclatt; // 4d coarse-coarse lattice + Coordinate c5simd, c5mpi; // 5D coarse simd/mpi + Coordinate cmsimd, cmmpi; // 6D coarse simd/mpi + // owned: + GridCartesian *Coarse5d; + GridCartesian *CoarseBatch; // 6D at the coarsening batch + GridCartesian *CoarseCoarse5d; + GridCartesian *CoarseCoarseBatch; + + MGCoarseGrids(GridCartesian *_FGrid, const MGSetupParams &P) + : FGrid(_FGrid) + { + Coordinate fdims = FGrid->FullDimensions(); // {Ls, x,y,z,t} + Coordinate fmpi = FGrid->_processors; + GRID_ASSERT( fdims.size() == 5 ); + Ls = fdims[0]; + batch = P.CoarsenBatch; + + clatt.resize(4); cclatt.resize(4); + for(int d=0;d<4;d++){ + GRID_ASSERT( fdims[d+1] % P.Block[d] == 0 ); + clatt[d] = fdims[d+1] / P.Block[d]; + } + for(int d=0;d<4;d++){ + GRID_ASSERT( clatt[d] % P.Block2[d] == 0 ); + cclatt[d] = clatt[d] / P.Block2[d]; + } + std::cout << GridLogMessage << "MGCoarseGrids: Block " << P.Block << " coarse lattice " << clatt << std::endl; + std::cout << GridLogMessage << "MGCoarseGrids: Block2 " << P.Block2 << " coarse-coarse lattice " << cclatt << std::endl; + + Coordinate c5latt({1,clatt[0],clatt[1],clatt[2],clatt[3]}); + c5simd = Coordinate({1,1,1,1,1}); + c5mpi = Coordinate({1,fmpi[1],fmpi[2],fmpi[3],fmpi[4]}); + Coarse5d = new GridCartesian(c5latt,c5simd,c5mpi); + + cmsimd = Coordinate({1,1,1,1,1,1}); + cmmpi = Coordinate({1,1,fmpi[1],fmpi[2],fmpi[3],fmpi[4]}); + Coordinate cblatt({batch,1,clatt[0],clatt[1],clatt[2],clatt[3]}); + CoarseBatch = new GridCartesian(cblatt,cmsimd,cmmpi); + + Coordinate cc5latt({1,cclatt[0],cclatt[1],cclatt[2],cclatt[3]}); + CoarseCoarse5d = new GridCartesian(cc5latt,c5simd,c5mpi); + + Coordinate ccblatt({batch,1,cclatt[0],cclatt[1],cclatt[2],cclatt[3]}); + CoarseCoarseBatch = new GridCartesian(ccblatt,cmsimd,cmmpi); + } + ~MGCoarseGrids() + { + delete CoarseCoarseBatch; + delete CoarseCoarse5d; + delete CoarseBatch; + delete Coarse5d; + } +}; + +////////////////////////////////////////////////////////////////////// +// ALL the coarsening information for one gauge configuration. +// +// Owns: the RAW near-null bases (fine and coarse -- retained so the +// coarse operators can be REBUILT on a changed gauge field with a +// fixed basis, the cheap HMC maintenance step; DiscardBasis() frees +// them for valence use), the two coarsened operators, the two block +// projectors, the dense bottom inverse (raw owning pointer: null in +// the constructor, allocated by BuildDenseBottom, deleted here), and +// the transient per-Nrhs solve grids (created by SetNrhs, which owns +// the ReleaseGrid -> delete -> new -> SetGrid ordering in ONE place). +// +// Borrows: the grid bundle. Declare MGCoarseGrids first. +////////////////////////////////////////////////////////////////////// +template +class PVdagMMultiGridCoarsening { +public: + typedef Lattice FineField; + typedef MultiGeneralCoarsenedOperatorV2 CoarseOperator; + typedef typename CoarseOperator::CoarseVector CoarseVector; + typedef typename CoarseVector::vector_object CoarseSiteObj; + typedef iScalar CComplex2; // coarsening deepens the nest by one iScalar + typedef MultiGeneralCoarsenedOperatorV2 CoarseCoarseOperator; + typedef typename CoarseCoarseOperator::CoarseVector CoarseCoarseVector; + typedef DenseCoarseMatrix DenseBottom; + + MGCoarseGrids &Grids; // borrowed + MGSetupParams Params; + NextToNearestStencilGeometry5D geom; + NextToNearestStencilGeometry5D geom2; + CoarseOperator CoarseOpPV; + CoarseCoarseOperator CoarseOpL2; + MultiRHSBlockProject MrhsProjector; + MultiRHSBlockProject MrhsProjectorL2; + DenseBottom *DenseCC; + std::vector rawNull; // RAW fine near-null basis + std::vector rawPsi; // RAW coarse near-null basis + GridCartesian *CMrhs; // transient solve grids, owned + GridCartesian *CCMrhs; + int nrhs; + + PVdagMMultiGridCoarsening(MGCoarseGrids &_Grids, const MGSetupParams &P) + : Grids(_Grids), + Params(P), + geom (_Grids.Coarse5d), + geom2(_Grids.CoarseCoarse5d), + CoarseOpPV(geom ,_Grids.Coarse5d), + CoarseOpL2(geom2,_Grids.CoarseCoarse5d), + DenseCC(nullptr), + CMrhs(nullptr), + CCMrhs(nullptr), + nrhs(-1) + {}; + + ~PVdagMMultiGridCoarsening() + { + if ( DenseCC ) delete DenseCC; + CoarseOpPV.ReleaseGrid(); // borrowers let go before their grids die + CoarseOpL2.ReleaseGrid(); + if ( CMrhs ) delete CMrhs; + if ( CCMrhs ) delete CCMrhs; + } + + //////////////////////////////////////////////////////////////////// + // The RAW fine basis: load from Params.SubspaceFile if it exists, + // else GCR inverse iteration from noise (and save if a file name + // was given). The Aggregation is scaffolding for CreateSubspaceGCR + // only -- that runs entirely on the fine grid, so the coarse grid + // it holds is never dereferenced. + //////////////////////////////////////////////////////////////////// + void GetSubspace(GridParallelRNG &RNG, LinearOperatorBase &FineOp) + { + uint64_t file_exists=0; + if ( Params.SubspaceFile.length() ) { + if ( Grids.FGrid->IsBoss() ){ std::ifstream f(Params.SubspaceFile); file_exists=f.good()?1:0; } + Grids.FGrid->GlobalSum(file_exists); + } + rawNull.clear(); + rawNull.reserve(nbasis); + for(int k=0;k Agg(Grids.Coarse5d,Grids.FGrid,0); + Agg.CreateSubspaceGCR(RNG,FineOp,nbasis); + for(int k=0;k + void Coarsen(PVdagMOp &FineOp) + { + GRID_ASSERT( rawNull.size() == nbasis ); + + // L1: working copy of the raw fine basis, orthonormalised in place + std::vector sub(rawNull); + + CoarseOpPV.SetGrid(Grids.CoarseBatch); + + std::cout << GridLogMessage << "PVdagMMultiGridCoarsening: L1 CoarsenOperator, batch " + << Grids.batch << std::endl; + + FineOp.SloppyComms(Params.FineSloppyComms); + CoarseOpPV.CoarsenOperator(FineOp,sub,Grids.Coarse5d,Grids.batch); + FineOp.SloppyComms(0); + + // Transfer operators from the orthonormalised basis; then the + // Galerkin images of the RAW basis define the L2 null space. + MrhsProjector.Allocate(nbasis,Grids.FGrid,Grids.Coarse5d); + MrhsProjector.ImportBasis(sub); + sub.clear(); sub.shrink_to_fit(); + + std::vector psi(nbasis,Grids.Coarse5d); + MrhsProjector.blockProject(rawNull,psi); + GramGuard("psi_coarse",psi,Grids.Coarse5d); + + rawPsi.clear(); + rawPsi.reserve(nbasis); + for(int k=0;k LinOpCoarse(CoarseOpPV); + std::cout << GridLogMessage << "PVdagMMultiGridCoarsening: L2 CoarsenOperator, batch " + << Grids.batch << std::endl; + CoarseOpL2.CoarsenOperator(LinOpCoarse,Grids.CoarseBatch,psi,Grids.CoarseCoarse5d); + + MrhsProjectorL2.Allocate(nbasis,Grids.Coarse5d,Grids.CoarseCoarse5d); + MrhsProjectorL2.ImportBasis(psi); // now block orthonormal + { + std::vector psi_cc(nbasis,Grids.CoarseCoarse5d); + MrhsProjectorL2.blockProject(rawPsi,psi_cc); // RAW vectors in + GramGuard("psi_cc",psi_cc,Grids.CoarseCoarse5d); + } + nrhs = -1; // both operators left at the batch grids + } + + //////////////////////////////////////////////////////////////////// + // Dense bottom: import the L2 matrix elements, invert (2D + // block-cyclic, fp64), then the check Import cannot run itself -- + // the mrhs operator applies only on a D+1 grid, so drive it at + // Nrhs 1 through a slice: ||A Ainv x - x||/||x||. + //////////////////////////////////////////////////////////////////// + void BuildDenseBottom(void) + { + if ( DenseCC ) delete DenseCC; + std::cout << GridLogMessage << "PVdagMMultiGridCoarsening: L3 dense bottom import" << std::endl; + DenseCC = new DenseBottom(Grids.CoarseCoarse5d); + DenseCC->Import(CoarseOpL2); + + Coordinate cc1latt({1,1,Grids.cclatt[0],Grids.cclatt[1],Grids.cclatt[2],Grids.cclatt[3]}); + GridCartesian *CoarseCoarseOne = new GridCartesian(cc1latt,Grids.cmsimd,Grids.cmmpi); + CoarseOpL2.SetGrid(CoarseCoarseOne); + + CoarseCoarseVector x(Grids.CoarseCoarse5d),y(Grids.CoarseCoarse5d),z(Grids.CoarseCoarse5d); + GridParallelRNG dRNG(Grids.CoarseCoarse5d); dRNG.SeedFixedIntegers(std::vector({11,12,13,14})); + random(dRNG,x); + + (*DenseCC)(x,y); // y = Ainv x + + CoarseCoarseVector y1(CoarseCoarseOne),z1(CoarseCoarseOne); + InsertSliceFast(y,y1,0,0); + CoarseOpL2.M(y1,z1); // z = A y + ExtractSliceFast(z,z1,0,0); + + z = z - x; + RealD rel = std::sqrt(norm2(z)/norm2(x)); + std::cout << GridLogMessage << "PVdagMMultiGridCoarsening: L3 dense ||A Ainv x - x||/||x|| = " + << rel << std::endl; + GRID_ASSERT( rel < 1.0e-2 ); + + CoarseOpL2.ReleaseGrid(); // let go before the grid dies + delete CoarseCoarseOne; + nrhs = -1; + } + + //////////////////////////////////////////////////////////////////// + // Retarget both operators to a solve Nrhs. The matrix elements are + // Nrhs independent and survive; only the D+1 grids and BLAS buffers + // change. The ordering discipline (borrowers release BEFORE their + // grids are destroyed) lives here and nowhere else. + //////////////////////////////////////////////////////////////////// + void SetNrhs(int nr) + { + if ( nr == nrhs ) return; + CoarseOpPV.ReleaseGrid(); + CoarseOpL2.ReleaseGrid(); + if ( CMrhs ) { delete CMrhs; CMrhs =nullptr; } + if ( CCMrhs ) { delete CCMrhs; CCMrhs=nullptr; } + Coordinate cml ({nr,1,Grids.clatt[0], Grids.clatt[1], Grids.clatt[2], Grids.clatt[3]}); + Coordinate ccml({nr,1,Grids.cclatt[0],Grids.cclatt[1],Grids.cclatt[2],Grids.cclatt[3]}); + CMrhs = new GridCartesian(cml, Grids.cmsimd,Grids.cmmpi); + CCMrhs = new GridCartesian(ccml,Grids.cmsimd,Grids.cmmpi); + CoarseOpPV.SetGrid(CMrhs); + CoarseOpL2.SetGrid(CCMrhs); + nrhs = nr; + std::cout << GridLogMessage << "PVdagMMultiGridCoarsening: operators at Nrhs " << nr << std::endl; + } + + //////////////////////////////////////////////////////////////////// + // Free the retained RAW bases (valence use: coarsening is final for + // this configuration and the fine basis is ~GB-scale host memory). + //////////////////////////////////////////////////////////////////// + void DiscardBasis(void) + { + rawNull.clear(); rawNull.shrink_to_fit(); + rawPsi.clear(); rawPsi.shrink_to_fit(); + } +}; + +////////////////////////////////////////////////////////////////////// +// The solve chain on a borrowed coarsening: adaptive PGCR smoothers +// at both levels (the one correct route for this non-Hermitian +// chain), dense bottom, mrhs V-cycle, outer shared-coefficient mrhs +// PGCR. Construction retargets the coarsening to this Nrhs. +// +// Halo policy: the whole V-cycle is preconditioner and runs with +// sloppy halos when FineSloppyComms is set; the outer Krylov is EXACT +// (its applications define what "converged" means), asserted by the +// FINAL true-residual report in Solve. +////////////////////////////////////////////////////////////////////// +template +class PVdagMMultiGridSolver { +public: + typedef typename Coarsening::FineField FineField; + typedef typename Coarsening::CoarseOperator CoarseOperator; + typedef typename Coarsening::CoarseVector CoarseVector; + typedef typename Coarsening::CoarseCoarseOperator CoarseCoarseOperator; + typedef typename Coarsening::CoarseCoarseVector CoarseCoarseVector; + typedef typename Coarsening::DenseBottom DenseBottom; + typedef PrecGeneralisedConjugateResidualNonHermitian FineSmoother_t; + typedef PrecGeneralisedConjugateResidualNonHermitian CoarseKrylov_t; + + Coarsening &C; + PVdagMMultiGridParams Params; + int nrhs; + int _regrid; // FIRST member-like init: SetNrhs before members capture grids + + PVdagMLinearOperator PVdagM; + ShiftedPVdagMLinearOperator ShiftedPVdagM; + TrivialPrecon simple_fine; + TrivialPrecon simpleC; + NonHermitianLinearOperator LinOpC; + ShiftedLinearOperator ShiftedC; + MrhsDenseCCSolve ccSolve; + CoarseKrylov_t CoarseSmootherGCR; + MrhsCoarseThreeLevelPrec L2to3Precon; + CoarseKrylov_t L2PGCR; + FineSmoother_t SmootherGCR; + MrhsTwoLevelMG ThreeLevelPrecon; + MrhsPGCRNonHermitian L1PGCR; + + PVdagMMultiGridSolver(Matrix &Ddwf, Matrix &Dpv, + Coarsening &_C, const PVdagMMultiGridParams &P, int nr) + : C(_C), Params(P), nrhs(nr), + _regrid((C.SetNrhs(nr),0)), // members below capture C.CMrhs/C.CCMrhs + PVdagM(Ddwf,Dpv), + ShiftedPVdagM(P.FineSmoother.Shift,Ddwf,Dpv), + LinOpC(C.CoarseOpPV), + ShiftedC(P.CoarseSmoother.Shift,LinOpC), + ccSolve(*C.DenseCC,nr), + // Coarse smoother: one preconditioned restart of nstep GCR steps. + CoarseSmootherGCR(0.01,1,ShiftedC,simpleC,P.CoarseSmoother.Mmax,P.CoarseSmoother.Nstep), + L2to3Precon(LinOpC,CoarseSmootherGCR,C.MrhsProjectorL2,ccSolve, + C.Grids.Coarse5d,C.Grids.CoarseCoarse5d,C.CCMrhs,nr), + // Coarse Krylov: Order/16 restarts of 16 steps. + L2PGCR(P.CoarseSolver.Tol,P.CoarseSolver.Order/16,LinOpC,L2to3Precon,P.CoarseSolver.Mmax,16), + // Fine smoother: one restart of nstep GCR steps, tolerance 0 = fixed work. + SmootherGCR(0.0,1,ShiftedPVdagM,simple_fine,P.FineSmoother.Mmax,P.FineSmoother.Nstep), + ThreeLevelPrecon(PVdagM,SmootherGCR,C.MrhsProjector,L2PGCR,C.Grids.Coarse5d,C.CMrhs), + L1PGCR(P.Outer.Tol,P.Outer.MaxIterations,PVdagM,ThreeLevelPrecon,P.Outer.Mmax,P.Outer.Nstep) + { + GRID_ASSERT( C.DenseCC != nullptr ); + + CoarseSmootherGCR.Level(2); + CoarseSmootherGCR.Name("Csmoother"); + CoarseSmootherGCR.SetZeroGuess(1); + + L2PGCR.Level(2); + L2PGCR.Name("Couter"); + L2PGCR.SetZeroGuess(1); + + SmootherGCR.Level(1); + SmootherGCR.Name("Fsmoother"); + SmootherGCR.SetZeroGuess(1); + + L1PGCR.Level(1); + L1PGCR.Name("Fouter"); + L1PGCR.SetZeroGuess(1); + + // The V-cycle is preconditioner: sloppy halos inside, exact restored on exit. + ThreeLevelPrecon.SetSloppy = [this](int s){ PVdagM.SloppyComms(s); }; + ThreeLevelPrecon.SloppyComms = Params.Setup.FineSloppyComms; + PVdagM.SloppyComms(0); + std::cout << GridLogMessage << "PVdagMMultiGridSolver: Nrhs " << nr + << ", fine halo policy: preconditioner+coarsening " + << (Params.Setup.FineSloppyComms ? "SLOPPY (fp32 wire)" : "exact") + << ", outer Krylov EXACT" << std::endl; + } + + void Solve(std::vector &src, std::vector &sol) + { + GRID_ASSERT( (int)src.size() == nrhs ); + // Start the solve from a clean device: evict setup-era Lattice copies + // and release the allocation caches' held blocks, so the LRU cap + // applies to the solve's own working set (OOM policy, Frontier). + MemoryManager::EvictAll(); + MemoryManager::DropCache(); + + GridStopWatch w; w.Start(); + L1PGCR(src,sol); + w.Stop(); + std::cout << GridLogMessage << "PVdagMMultiGridSolver: Nrhs "< + + 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. + + See the full license in the file "LICENSE" in the top level distribution + directory +*************************************************************************************/ +/* END LEGAL */ +#pragma once + +NAMESPACE_BEGIN(Grid); + +////////////////////////////////////////////////////////////////////////////////////// +// Serialisable parameters for the mrhs PVdagM multigrid solver chain. +// +// Named for the multigrid FORM it configures, as WilsonMGParams is for the +// Wilson MG tests; an MrhsHDCGParams will compose the same generic MG* +// sub-structs for the Hermitian chain. Read from XML (or JSON -- the +// serialisation macros give both) via ReadPVdagMMultiGridParams below and +// printed at startup by the macro's operator<<, so every log identifies its +// own run. The constructor defaults ARE the banked optimum; an +// unconfigured run reproduces the best recorded point. Update them when a +// better point is banked, and date the change. +// +// Current optimum: 2026-08-24, slurm-5335492 F4, 48^3x96 Ls=24 on 288 GCDs, +// 17.2 s/RHS at Nrhs=4, 32.2 s at Nrhs=1 (exact-halo FINAL ~1e-8). +// Smoother mmax == nstep (full GCR history). +// +// There are NO environment-variable controls anywhere in this subsystem +// (ruling 2026-09-05): parameters live here, library-internal constants are +// hard defaults in their classes, verbosity is the --log channels. +// Research instruments (GCR coefficient recording/replay, Chebyshev and +// stationary smoother variants, power iteration) are deliberately NOT part +// of this interface: they remain programmatic, for algorithmic studies, +// not consumer knobs. +////////////////////////////////////////////////////////////////////////////////////// + +// The smoother is the adaptive shifted PGCR -- the one correct route for +// this non-Hermitian chain (stationary replay and Chebyshev were explored +// and did not win here; Chebyshev remains the HERMITIAN chain's smoother). +struct MGSmootherParams : Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(MGSmootherParams, + RealD, Shift, + int, Nstep, + int, Mmax); + MGSmootherParams(RealD shift=0.1, int nstep=6, int mmax=6) + : Shift(shift), Nstep(nstep), Mmax(mmax) {}; +}; + +struct MGCoarseSolverParams : Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(MGCoarseSolverParams, + RealD, Tol, + int, Order, + int, Mmax); + MGCoarseSolverParams() : Tol(0.05), Order(200), Mmax(16) {}; +}; + +struct MGOuterParams : Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(MGOuterParams, + RealD, Tol, + int, MaxIterations, + int, Mmax, + int, Nstep); + MGOuterParams() : Tol(1.0e-8), MaxIterations(1000), Mmax(4), Nstep(8) {}; +}; + +struct MGDenseParams : Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(MGDenseParams, + int, LeafSpan); // big-leaf span in blocks (banked 9); see BlockCyclicSchurInverse + MGDenseParams() : LeafSpan(9) {}; +}; + +struct MGSetupParams : Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(MGSetupParams, + std::vector, Block, // fine -> coarse blocking + std::vector, Block2, // coarse -> coarse-coarse blocking + int, CoarsenBatch, + std::string, SubspaceFile, // scidac; empty = create from noise, no I/O + int, FineSloppyComms);// fp32 wire INSIDE the preconditioner only + MGSetupParams() + : Block({2,2,3,3}), Block2({4,4,2,4}), CoarsenBatch(9), + SubspaceFile(""), FineSloppyComms(1) {}; +}; + +struct PVdagMMultiGridParams : Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(PVdagMMultiGridParams, + MGSetupParams, Setup, + MGSmootherParams, FineSmoother, + MGSmootherParams, CoarseSmoother, + MGCoarseSolverParams, CoarseSolver, + MGOuterParams, Outer, + MGDenseParams, Dense); + PVdagMMultiGridParams() + : FineSmoother (0.1, 6, 6), + CoarseSmoother(0.1, 2, 2) {}; +}; + +inline void CheckValidity(const PVdagMMultiGridParams &P) +{ + GRID_ASSERT( P.Setup.Block.size() == 4 ); + GRID_ASSERT( P.Setup.Block2.size() == 4 ); + GRID_ASSERT( P.Setup.CoarsenBatch >= 1 ); + GRID_ASSERT( P.FineSmoother.Nstep > 0 ); + GRID_ASSERT( P.CoarseSmoother.Nstep > 0 ); + GRID_ASSERT( P.CoarseSolver.Tol > 0.0 ); + GRID_ASSERT( P.Outer.Tol > 0.0 ); + GRID_ASSERT( P.Dense.LeafSpan >= 1 ); +} + +////////////////////////////////////////////////////////////////////////////////////// +// Read from an XML file; on a MISSING file the boss writes .templ from +// the defaults and returns false so the caller can Grid_finalize() and exit +// -- the produced template is the documentation of every parameter. An +// empty filename means "defaults": no I/O, returns true. Always prints the +// effective parameters, so the log describes the run. +// +// Drivers fetch the filename themselves, e.g. +// std::string pfile(""); +// if( GridCmdOptionExists(argv,argv+argc,"--multigrid-params") ) +// pfile = GridCmdOptionPayload(argv,argv+argc,"--multigrid-params"); +// if( !ReadPVdagMMultiGridParams(Params, pfile) ) { Grid_finalize(); return 0; } +////////////////////////////////////////////////////////////////////////////////////// +inline bool ReadPVdagMMultiGridParams(PVdagMMultiGridParams &P, const std::string &xmlfile) +{ + if ( xmlfile.length() ) { + bool good; + { std::ifstream f(xmlfile); good = f.good(); } + if ( !good ) { + std::cout << GridLogMessage << "PVdagMMultiGridParams: " << xmlfile << " does not exist" << std::endl; + if ( GlobalSharedMemory::WorldRank == 0 ) { + XmlWriter WR(xmlfile+".templ"); + write(WR, "PVdagMMultiGridParams", P); + std::cout << GridLogMessage << "PVdagMMultiGridParams: template written to " + << xmlfile << ".templ" << std::endl; + } + return false; + } + XmlReader RD(xmlfile); + read(RD, "PVdagMMultiGridParams", P); + } + CheckValidity(P); + std::cout << GridLogMessage << "PVdagMMultiGridParams (" + << (xmlfile.length() ? xmlfile : std::string("defaults")) << "):" << std::endl; + std::cout << P << std::endl; + return true; +} + +NAMESPACE_END(Grid); diff --git a/Grid/algorithms/multigrid/PVdagMOperators.h b/Grid/algorithms/multigrid/PVdagMOperators.h new file mode 100644 index 000000000..6e2bef4ed --- /dev/null +++ b/Grid/algorithms/multigrid/PVdagMOperators.h @@ -0,0 +1,97 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./Grid/algorithms/multigrid/PVdagMOperators.h + + Copyright (C) 2026 + +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. + + See the full license in the file "LICENSE" in the top level distribution + directory +*************************************************************************************/ +/* END LEGAL */ +#pragma once + +NAMESPACE_BEGIN(Grid); + +////////////////////////////////////////////////////////////////////// +// A = PV^dag M (non-Hermitian). The operator the PVdagM multigrid +// coarsens and solves. +// +// SloppyComms delegates to BOTH underlying fermion operators, so one +// call controls every fine halo in the chain -- including any shifted +// wrapper built on the same two matrices. Halo-precision POLICY: +// reduced-precision (fp32 wire) halos belong in the PRECONDITIONER -- +// the smoother, the V-cycle's own residuals, the coarsening -- and +// NEVER in the outer Krylov, whose applications define what +// "converged" means. The operators default to EXACT; the V-cycle +// turns sloppiness on for its own scope only (MrhsTwoLevelMG). +////////////////////////////////////////////////////////////////////// +template +class PVdagMLinearOperator : public LinearOperatorBase { + Matrix &_Mat; Matrix &_PV; +public: + PVdagMLinearOperator(Matrix &Mat,Matrix &PV): _Mat(Mat),_PV(PV) {}; + void SloppyComms(int sloppy) { _Mat.SloppyComms(sloppy); _PV.SloppyComms(sloppy); } + void OpDiag (const Field &in, Field &out) { GRID_ASSERT(0); } + void OpDir (const Field &in, Field &out,int dir,int disp) { GRID_ASSERT(0); } + void OpDirAll (const Field &in, std::vector &out){ GRID_ASSERT(0); }; + void Op (const Field &in, Field &out){ Field tmp(in.Grid()); _Mat.M(in,tmp); _PV.Mdag(tmp,out); } + void AdjOp (const Field &in, Field &out){ Field tmp(in.Grid()); _PV.M(in,tmp); _Mat.Mdag(tmp,out); } + void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ HermOp(in,out); ComplexD d=innerProduct(in,out); n1=real(d); n2=norm2(out); } + void HermOp(const Field &in, Field &out){ Field tmp(in.Grid()); Op(in,tmp); AdjOp(tmp,out); } +}; + +////////////////////////////////////////////////////////////////////// +// A + shift, built directly on the two matrices (fine level: avoids +// an extra fine-field pass through a generic shifted wrapper). +////////////////////////////////////////////////////////////////////// +template +class ShiftedPVdagMLinearOperator : public LinearOperatorBase { + Matrix &_Mat; Matrix &_PV; +public: + RealD shift; + ShiftedPVdagMLinearOperator(RealD _shift,Matrix &Mat,Matrix &PV): shift(_shift),_Mat(Mat),_PV(PV){}; + void OpDiag (const Field &in, Field &out) { GRID_ASSERT(0); } + void OpDir (const Field &in, Field &out,int dir,int disp) { GRID_ASSERT(0); } + void OpDirAll (const Field &in, std::vector &out){ GRID_ASSERT(0); }; + void Op (const Field &in, Field &out){ Field tmp(in.Grid()); _Mat.M(in,tmp); _PV.Mdag(tmp,out); out = out + shift*in; } + // BUG LEDGER 2026-09-06: the example-local original read + // _PV.M(tmp,out); _Mat.Mdag(in,tmp); + // consuming tmp before writing it -- latent, never exercised (nothing in + // the chain calls AdjOp on the shifted fine operator). Corrected here to + // the adjoint of Op, matching PVdagMLinearOperator::AdjOp plus the shift. + void AdjOp (const Field &in, Field &out){ Field tmp(in.Grid()); _PV.M(in,tmp); _Mat.Mdag(tmp,out); out = out + shift*in; } + void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ GRID_ASSERT(0); } + void HermOp(const Field &in, Field &out){ Field tmp(in.Grid()); Op(in,tmp); AdjOp(tmp,out); } +}; + +////////////////////////////////////////////////////////////////////// +// Op + shift for any operator (coarse levels). +////////////////////////////////////////////////////////////////////// +template +class ShiftedLinearOperator : public LinearOperatorBase { + LinearOperatorBase &_Op; RealD shift; +public: + ShiftedLinearOperator(RealD _shift, LinearOperatorBase &Op) : _Op(Op), shift(_shift) {} + void OpDiag (const Field &in, Field &out) { GRID_ASSERT(0); } + void OpDir (const Field &in, Field &out,int dir,int disp) { GRID_ASSERT(0); } + void OpDirAll (const Field &in, std::vector &out) { GRID_ASSERT(0); } + void Op (const Field &in, Field &out) { _Op.Op(in,out); out = out + shift*in; } + void AdjOp (const Field &in, Field &out) { _Op.AdjOp(in,out); out = out + shift*in; } + void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ GRID_ASSERT(0); } + void HermOp (const Field &in, Field &out) { Field tmp(in.Grid()); Op(in,tmp); AdjOp(tmp,out); } +}; + +// The non-Hermitian spectral-edge diagnostic that used to live here as +// PowerIteration() now lives beside its Hermitian sibling as +// Grid::NonHermitianPowerMethod (Grid/algorithms/iterative/PowerMethod.h). + +NAMESPACE_END(Grid);