diff --git a/Grid/algorithms/multigrid/HDCGMultiGrid.h b/Grid/algorithms/multigrid/HDCGMultiGrid.h new file mode 100644 index 000000000..065d4bd9f --- /dev/null +++ b/Grid/algorithms/multigrid/HDCGMultiGrid.h @@ -0,0 +1,632 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./Grid/algorithms/multigrid/HDCGMultiGrid.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 +#include +#include +#include + +NAMESPACE_BEGIN(Grid); + +////////////////////////////////////////////////////////////////////////////////////// +// The two-level mrhs HDCG (arXiv:1402.2585, mrhs form arXiv:2409.03904) as +// objects, the shape of PVdagMMultiGrid.h: +// +// MGCoarseGrids / MGFineGridsF shared with the PVdagM chain +// HDCGCoarsening everything that is a function of the gauge +// field: the raw near-null basis (and its +// refinement), the transfer operators, the +// Galerkin coarse operator (V2, Hermitian), +// the coarse eigenvectors for deflation +// HDCGSolver the solve chain on a borrowed coarsening: +// smoother, deflated coarse solve, ADEF-2 +// preconditioner (fp64 or fp32 behind the +// seam), outer mrhs fPcg / BlockCGrQ +// +// The Hermitian operator is the Schur-preconditioned M^dag M on one +// checkerboard; the near-null basis and the fine fields live on the +// checkerboarded grid. The coarse space is unvectorised and ALWAYS the D+1 +// multiRHS field; Nrhs 1 is the extent-1 case. Grids outlive the +// coarsening outlives the solver. +////////////////////////////////////////////////////////////////////////////////////// +template +class HDCGCoarsening { +public: + typedef Lattice FineField; + typedef MultiGeneralCoarsenedOperatorV2 CoarseOperator; + typedef typename CoarseOperator::CoarseVector CoarseVector; + typedef typename GridTypeMapper::SinglePrecision FobjF; + typedef Lattice FineFieldF; + + MGCoarseGrids &Grids; // borrowed + MGFineGridsF &GridsF; // borrowed + GridBase *FrbGrid; // borrowed: the checkerboarded fine grid + int cb; + HDCGMultiGridParams Params; + NextToNextToNextToNearestStencilGeometry5D geom; // the coarse stencil of the Schur operator + CoarseOperator CoarseOp; + ////////////////////////////////////////////////////////////////////// + // ONE fine transfer operator. Its STORE follows the coarse sector -- + // the coarse space is what it feeds -- while its import and export + // accept EITHER fine precision, because they are already a layout + // transformation and a scalar conversion inside one costs nothing. + ////////////////////////////////////////////////////////////////////// + static const bool CoarseIsSingle = + ( sizeof(typename GridTypeMapper::scalar_type) == sizeof(ComplexF) ); + typedef typename std::conditional::type ProjectorField; + typedef MultiRHSBlockProject ProjectorL1_t; + + ProjectorL1_t MrhsProjectorL1; + std::vector rawNull; // RAW fine near-null basis + std::vector evec; // coarse eigenvectors, on Coarse5d + std::vector eval; + MultiRHSDeflation Deflator; // nev==0 until CoarseLanczos + GridCartesian *CMrhs; // transient solve grid, owned + int nrhs; + + HDCGCoarsening(MGCoarseGrids &_Grids, MGFineGridsF &_GridsF, GridBase *_FrbGrid, int _cb, + const HDCGMultiGridParams &P) + : Grids(_Grids), GridsF(_GridsF), FrbGrid(_FrbGrid), cb(_cb), Params(P), + geom(_Grids.Coarse5d), + CoarseOp(geom,_Grids.Coarse5d), + CMrhs(nullptr), + nrhs(-1) + { + Deflator.Deallocate(); + // The coarse stencil is a max-norm-1 box (Geometry.h) and the Schur + // operator M^dag M reaches four fine hops, so a block shorter than 4 + // leaves couplings outside the stencil: the Fourier extraction then + // folds them into the retained points and the Galerkin certificate + // reads it. Warn here; the certificate is the guard. + for(int d=0;d<4;d++){ + if ( Params.Setup.Block1[d] < 4 ) + std::cout << GridLogWarning << "HDCGCoarsening: Block1[" << d << "] = " << Params.Setup.Block1[d] + << " < 4: the Schur operator's range exceeds the coarse stencil; " + << "expect a Galerkin defect" << std::endl; + } + }; + ~HDCGCoarsening() + { + CoarseOp.ReleaseGrid(); // borrowers let go before their grids die + if ( CMrhs ) delete CMrhs; + } + + int FileExists(const std::string &fname) + { + uint64_t exists=0; + if ( fname.length() ) { + if ( Grids.FGrid->IsBoss() ){ std::ifstream f(fname); exists=f.good()?1:0; } + Grids.FGrid->GlobalSum(exists); + } + return (int)exists; + } + FineField FineNoise(GridParallelRNG &RNG) + { + FineField v(FrbGrid); v.Checkerboard()=cb; gaussian(RNG,v); + return v; + } + // Chebyshev upper bound: the parameter, or 1.1 x a power-method estimate + RealD UpperBound(RealD param, LinearOperatorBase &HermOp, GridParallelRNG &RNG) + { + if ( param > 0.0 ) return param; + FineField src = FineNoise(RNG); + PowerMethod PM; + RealD lam = PM(HermOp,src); + std::cout << GridLogMessage << "HDCGCoarsening: fine power method " << lam + << ", Chebyshev upper bound " << 1.1*lam << std::endl; + return 1.1*lam; + } + + //////////////////////////////////////////////////////////////////// + // The RAW basis: a refined basis on disk wins, then a raw basis on + // disk, else built by the chosen method (and saved if named). CG + // refinement happens here; HDCG refinement needs a solver, see + // HDCGRefineSubspace below. The Aggregation is scaffolding for the + // Chebyshev filters only: they run entirely on the fine grid. + //////////////////////////////////////////////////////////////////// + void GetSubspace(GridParallelRNG &RNG, LinearOperatorBase &HermOp) + { + const HDCGSubspaceParams &S = Params.Subspace; + rawNull.clear(); + rawNull.reserve(nbasis); + for(int k=0;k Agg(Grids.Coarse5d,FrbGrid,cb); + if ( S.Method == HDCGSubspaceMethod::ChebyshevNew ) { + RealD hi = UpperBound(S.ChebyHi,HermOp,RNG); + std::cout << GridLogMessage << "HDCGCoarsening: CreateSubspaceChebyshevNew hi " << hi << std::endl; + Agg.CreateSubspaceChebyshevNew(RNG,HermOp,hi); + } else if ( S.Method == HDCGSubspaceMethod::Chebyshev ) { + RealD hi = UpperBound(S.ChebyHi,HermOp,RNG); + std::cout << GridLogMessage << "HDCGCoarsening: CreateSubspaceChebyshev hi " << hi + << " lo " << S.ChebyLo << " order " << S.ChebyOrder << std::endl; + Agg.CreateSubspaceChebyshev(RNG,HermOp,nbasis,hi,S.ChebyLo,S.ChebyOrder); + } else { + std::cout << GridLogMessage << "HDCGCoarsening: subspace from fine eigenvectors " << S.FineEvecFile + << " summed " << S.FineEvecSample << " per vector" << std::endl; + LoadFineEvecSum(Agg.subspace,S.FineEvecFile,S.FineEvecSample,S.FileControl); + } + for(int k=0;k Agg(Grids.Coarse5d,FrbGrid,cb); + for(int k=0;k &sub, const std::string &file, int sample, int control) + { +#ifdef HAVE_LIME + emptyUserRecord record; + FineFieldF tmp(GridsF.FrbGridF), sum(GridsF.FrbGridF); + tmp.Checkerboard()=cb; sum.Checkerboard()=cb; + ScidacReader RD; + RD.open(file); + for(int b=0;b &HermOp, TwoLevelCGmrhs &Refiner, int nr) + { + std::vector src(nr,FrbGrid), res(nr,FrbGrid); + FineField tmp(FrbGrid); tmp.Checkerboard()=cb; + for(int r=0;r " << norm2(tmp) << std::endl; + Refiner(src,res); + for(int r=0;r " << norm2(tmp) << std::endl; + } + } + + //////////////////////////////////////////////////////////////////// + // The Galerkin coarse operator and transfer operators from the RAW + // basis. CoarsenOperator block-orthonormalises its input in place, + // so a working copy is taken and the RAW vectors survive. The + // operator handed in must apply the Hermitian operator through Op + // (HermOpAdaptor). + // + // The RAW basis has done its work once this returns, so it is freed -- + // it is nbasis fine vectors, the single largest resident block of the + // setup. retain_basis keeps it, for the two callers that coarsen more + // than once: the HDCG subspace refinement below, and a fixed-basis + // rebuild on a changed gauge field (HMC). Setup.RetainSubspace keeps + // it for the whole life of the object. + //////////////////////////////////////////////////////////////////// + void Coarsen(LinearOperatorBase &HermOp, int retain_basis=0) + { + GRID_ASSERT( rawNull.size() == nbasis ); + std::vector sub(rawNull); + + CoarseOp.SetGrid(Grids.CoarseBatch); + std::cout << GridLogMessage << "HDCGCoarsening: CoarsenOperator, batch " << Grids.batch << std::endl; + // CoarsenOperator orthonormalises sub in place and imports it into OUR + // transfer operator, which it then uses. One basis store, imported once. + // The projector's grid carries the STORE's SIMD layout: the fp32 fine + // red-black grid when the coarse sector is fp32. + MrhsProjectorL1.Allocate(nbasis, + CoarseIsSingle ? (GridBase *)GridsF.FrbGridF : (GridBase *)FrbGrid, + Grids.Coarse5d); + CoarseOp.CoarsenOperator(HermOp,sub,Grids.Coarse5d,Grids.batch,MrhsProjectorL1); + if ( !retain_basis && !Params.Setup.RetainSubspace ) DiscardBasis(); + MrhsProjectorL1.ReleaseScratch(); + nrhs = -1; // the batch grid is set; SetNrhs retargets + ReportDeviceFootprint("after Coarsen"); + } + + //////////////////////////////////////////////////////////////////// + // Retarget the coarse operator to Nrhs: the D+1 grid is owned here + // and every borrower re-captures it after this call. + //////////////////////////////////////////////////////////////////// + void SetNrhs(int nr) + { + if ( nr == nrhs ) return; + CoarseOp.ReleaseGrid(); + if ( CMrhs ) { delete CMrhs; CMrhs=nullptr; } + Coordinate cml({nr,1,Grids.clatt[0],Grids.clatt[1],Grids.clatt[2],Grids.clatt[3]}); + CMrhs = new GridCartesian(cml,Grids.cmsimd,Grids.cmmpi); + CoarseOp.SetGrid(CMrhs); + nrhs = nr; + std::cout << GridLogMessage << "HDCGCoarsening: coarse operator at Nrhs " << nr << std::endl; + } + + //////////////////////////////////////////////////////////////////// + // Coarse eigenvectors for the deflated coarse solve. Block Lanczos + // runs at the solve's Nrhs on the D+1 operator; plain Lanczos at + // Nrhs 1 (the extent-1 D+1 grid), its vectors then sliced to the D + // grid so both feed the deflator the same way. Loaded from EvecFile + // when present, saved there when named. + //////////////////////////////////////////////////////////////////// + void CoarseLanczos(int nr) + { + const HDCGLanczosParams &L = Params.Lanczos; + evec.clear(); eval.clear(); + Deflator.Deallocate(); + if ( L.Method == HDCGLanczosMethod::None ) { + std::cout << GridLogMessage << "HDCGCoarsening: no coarse deflation" << std::endl; + return; + } + if ( FileExists(L.EvecFile) ) { + std::cout << GridLogMessage << "HDCGCoarsening: loading coarse eigenpairs " << L.EvecFile << std::endl; + loadEigenpairs(evec,eval,L.EvecFile,Grids.Coarse5d,Params.Subspace.FileControl); + } else { + int Nconv=0; + SetNrhs( L.Method == HDCGLanczosMethod::Block ? nr : 1 ); + HermitianLinearOperator HermOpC(CoarseOp); + GridParallelRNG cRNG(CMrhs); cRNG.SeedFixedIntegers(std::vector({31,32,33,34})); + RealD hi = L.ChebyHi; + if ( hi <= 0.0 ) { + CoarseVector s(CMrhs); random(cRNG,s); + PowerMethod PM; + hi = 1.1*PM(HermOpC,s); + std::cout << GridLogMessage << "HDCGCoarsening: coarse Chebyshev upper bound " << hi << std::endl; + } + Chebyshev Cheby(L.ChebyLo,hi,L.ChebyOrder); + + if ( L.Method == HDCGLanczosMethod::Block ) { + std::cout << GridLogMessage << "HDCGCoarsening: block Lanczos, Nrhs " << nr + << " Nstop " << L.Nstop << " Nk " << L.Nk << " Nm " << L.Nm << std::endl; + ImplicitlyRestartedBlockLanczosCoarse + IRL(HermOpC,Grids.Coarse5d,CMrhs,nr,Cheby,L.Nstop,1,nr,L.Nk,L.Nm,L.Tol,L.MaxIt); + evec.resize(L.Nm,Grids.Coarse5d); + eval.resize(L.Nm); + std::vector src(nr,Grids.Coarse5d); + GridParallelRNG dRNG(Grids.Coarse5d); dRNG.SeedFixedIntegers(std::vector({41,42,43,44})); + for(int r=0;r Op(HermOpC); + FunctionHermOp OpCheby(Cheby,HermOpC); + ImplicitlyRestartedLanczos IRL(OpCheby,Op,L.Nstop,L.Nk,L.Nm,L.Tol,L.MaxIt); + std::vector evec1(L.Nm,CMrhs); + std::vector eval1(L.Nm); + CoarseVector src(CMrhs); random(cRNG,src); + IRL.calc(eval1,evec1,src,Nconv); + evec.resize(Nconv,Grids.Coarse5d); + eval.resize(Nconv); + for(int e=0;e - | / ||, which CG relies on. + //////////////////////////////////////////////////////////////////// + RealD CertifyCoarsening(LinearOperatorBase &HermOp) + { + SetNrhs(1); + CoarseVector xc(CMrhs), yc(CMrhs), Acx(CMrhs), Acy(CMrhs), PtAPx(CMrhs); + GridParallelRNG cRNG(CMrhs); cRNG.SeedFixedIntegers(std::vector({21,22,23,24})); + random(cRNG,xc); + random(cRNG,yc); + + CoarseOp.M(xc,Acx); // A_c x + CoarseOp.M(yc,Acy); // A_c y + + std::vector Px(1,FrbGrid), APx(1,FrbGrid); + Px[0].Checkerboard()=cb; APx[0].Checkerboard()=cb; + MrhsProjectorL1.blockPromote(Px,xc); // P x + HermOp.HermOp(Px[0],APx[0]); // A P x + MrhsProjectorL1.blockProject(APx,PtAPx); // P^dag A P x + + CoarseVector d(CMrhs); d = Acx - PtAPx; + RealD rel = std::sqrt(norm2(d)/norm2(Acx)); + ComplexD xAy = TensorRemove(innerProduct(xc,Acy)); + ComplexD Axy = TensorRemove(innerProduct(Acx,yc)); + RealD herm = abs(xAy-Axy)/abs(xAy); + std::cout << GridLogMessage << "HDCGCoarsening: GALERKIN CERTIFICATE ||A_c x - P^dag A P x||/||A_c x|| = " << rel << std::endl; + std::cout << GridLogMessage << "HDCGCoarsening: HERMITICITY CERTIFICATE | - |/|| = " << herm << std::endl; + return rel; + } + + + ////////////////////////////////////////////////////////////////////// + // The device memory this coarsening holds that the memory manager + // CANNOT evict: BLAS stores, which are plain device allocations. + // Lattice fields are omitted on purpose -- they are evictable, so they + // cost eviction traffic, not failure. What is reported is what has to + // fit alongside the manager's cache, the comms buffers and the shm + // segment, which is the budget an out-of-memory run actually breaks. + ////////////////////////////////////////////////////////////////////// + void ReportDeviceFootprint(const std::string &when) + { + uint64_t p1 = MrhsProjectorL1.DeviceBytes(); + uint64_t o1 = CoarseOp.DeviceBytes(); + uint64_t df = Deflator.DeviceBytes(); + const double G = 1024.*1024.*1024.; + std::cout << GridLogMessage << "HDCGCoarsening (" << when << "): resident device memory/rank " + << (p1+o1+df)/G << " GB = transfer L1 " << p1/G + << " + operator L1 " << o1/G + << " + deflation " << df/G << " GB" << std::endl; + std::cout << GridLogMessage << "HDCGCoarsening (" << when << "): not evictable; must fit alongside the " + << MemoryManager::DeviceMaxBytes/G + << " GB manager cache, the comms buffers and the shm segment" << std::endl; + } + + void DiscardBasis(void) + { + rawNull.clear(); rawNull.shrink_to_fit(); + } +}; + +////////////////////////////////////////////////////////////////////// +// The solve chain on a borrowed coarsening. Both fine precisions are +// built; the outer picks the fp32 preconditioner through the seam when +// Setup.FinePrecision is fp32. OperatorShift is 0 in production; the +// HDCG-refinement solver runs on the shifted operator. +////////////////////////////////////////////////////////////////////// +template +class HDCGSolver { +public: + typedef typename Coarsening::FineField FineField; + typedef typename Coarsening::FineFieldF FineFieldF; + typedef typename Coarsening::CoarseOperator CoarseOperator; + typedef typename Coarsening::CoarseVector CoarseVector; + + Coarsening &C; + HDCGMultiGridParams Params; + int nrhs; + int _regrid; // FIRST member-like init: SetNrhs before members capture grids + RealD OperatorShift; + + // fp64 fine chain + SchurDiagMooeeOperator HermOpEO; + ShiftedHermOpLinearOperator HermOpShifted; + HermOpAdaptor FineOp; // Op = HermOp: the operator of the solve + ShiftedHermOpLinearOperator SmootherOp; + CGSmoother SmootherCG; + ChebyshevSmoother SmootherCheby; + // coarse + HermitianLinearOperator LinOpC; + ConjugateGradient CoarseCG; + DoNothingGuesser NoGuess; + HPDSolver CoarseDeflatedCG; + ChebyshevInverter CoarseCheby; + MrhsADEF2Preconditioner ADEF2; + // fp32 fine chain, behind the seam; the coarse solve is shared + SchurDiagMooeeOperator HermOpEOF; + ShiftedHermOpLinearOperator HermOpShiftedF; + HermOpAdaptor FineOpF; + ShiftedHermOpLinearOperator SmootherOpF; + CGSmoother SmootherCGF; + ChebyshevSmoother SmootherChebyF; + MrhsADEF2Preconditioner ADEF2F; + MrhsMixedPrecPreconditioner PreconSeam; + TwoLevelCGmrhs Outer; + + LinearFunction &PickSmoother(void) { + return (Params.Smoother.Type==HDCGSmootherType::CG) + ? static_cast&>(SmootherCG) + : static_cast&>(SmootherCheby); + } + LinearFunction &PickSmootherF(void) { + return (Params.Smoother.Type==HDCGSmootherType::CG) + ? static_cast&>(SmootherCGF) + : static_cast&>(SmootherChebyF); + } + LinearFunction &PickCoarseSolve(void) { + return (Params.CoarseSolver.Method==HDCGCoarseMethod::DeflatedCG) + ? static_cast&>(CoarseDeflatedCG) + : static_cast&>(CoarseCheby); + } + template + static LinearOperatorBase &Unshifted(RealD shift, Op &op, ShiftedOp &shifted) { + return (shift==0.0) ? static_cast&>(op) : static_cast&>(shifted); + } + + HDCGSolver(Matrix &Ddwf, MatrixF &DdwfF, Coarsening &_C, const HDCGMultiGridParams &P, int nr, RealD shift=0.0) + : C(_C), Params(P), nrhs(nr), + _regrid((C.SetNrhs(nr),0)), + OperatorShift(shift), + HermOpEO(Ddwf), + HermOpShifted(HermOpEO,shift), + FineOp(Unshifted(shift,HermOpEO,HermOpShifted)), + SmootherOp(FineOp,P.Smoother.Shift), + SmootherCG(P.Smoother.Order,SmootherOp), + SmootherCheby(P.Smoother.ChebyLo,P.Smoother.ChebyHi,P.Smoother.Order,FineOp), + LinOpC(C.CoarseOp), + CoarseCG(P.CoarseSolver.Tol,P.CoarseSolver.MaxIt,false), + CoarseDeflatedCG(LinOpC,CoarseCG,NoGuess), + CoarseCheby(P.CoarseSolver.ChebyLo,P.CoarseSolver.ChebyHi,P.CoarseSolver.ChebyOrder,LinOpC), + ADEF2(FineOp,PickSmoother(),PickCoarseSolve(),PickCoarseSolve(),C.MrhsProjectorL1,C.Deflator,C.CMrhs,C.FrbGrid), + HermOpEOF(DdwfF), + HermOpShiftedF(HermOpEOF,shift), + FineOpF(Unshifted(shift,HermOpEOF,HermOpShiftedF)), + SmootherOpF(FineOpF,P.Smoother.Shift), + SmootherCGF(P.Smoother.Order,SmootherOpF), + SmootherChebyF(P.Smoother.ChebyLo,P.Smoother.ChebyHi,P.Smoother.Order,FineOpF), + ADEF2F(FineOpF,PickSmootherF(),PickCoarseSolve(),PickCoarseSolve(),C.MrhsProjectorL1,C.Deflator,C.CMrhs,C.GridsF.FrbGridF), + PreconSeam(ADEF2F,C.FrbGrid,C.GridsF.FrbGridF,nr), + Outer(P.Outer.Tol,P.Outer.MaxIt,FineOp, + (P.Setup.FinePrecision==MGPrecision::fp32) + ? static_cast&>(PreconSeam) + : static_cast&>(ADEF2), + C.FrbGrid, + (P.Outer.Algorithm==HDCGOuterAlgorithm::PrecBlockCGrQ) ? MrhsCGAlgorithm::PrecBlockCGrQ : MrhsCGAlgorithm::fPcg) + { + // BlockCGrQ assumes a stationary preconditioner: the fixed-iteration CG + // smoother is fixed WORK but a nonlinear map, and a coarse CG solved to a + // tolerance is not stationary either. Chebyshev on both is. + if ( Params.Outer.Algorithm==HDCGOuterAlgorithm::PrecBlockCGrQ && + ( Params.CoarseSolver.Method!=HDCGCoarseMethod::Chebyshev || + Params.Smoother.Type !=HDCGSmootherType::Chebyshev ) ) + std::cout << GridLogWarning << "HDCGSolver: PrecBlockCGrQ with a non-stationary preconditioner " + << "(needs Smoother.Type Chebyshev AND CoarseSolver.Method Chebyshev); expect stagnation" + << std::endl; + std::cout << GridLogMessage << "HDCGSolver: Nrhs " << nr + << ", operator shift " << shift + << ", smoother " << (Params.Smoother.Type==HDCGSmootherType::CG ? "CG" : "Chebyshev") + << " order " << Params.Smoother.Order + << ", coarse " << (Params.CoarseSolver.Method==HDCGCoarseMethod::DeflatedCG ? "deflated CG" : "Chebyshev") + << " (deflation vectors " << C.Deflator.nev << ")" + << ", preconditioner fine level " + << (Params.Setup.FinePrecision==MGPrecision::fp32 ? "fp32 (seam at the outer CG)" : "fp64") + << ", outer " << (Params.Outer.Algorithm==HDCGOuterAlgorithm::PrecBlockCGrQ ? "PrecBlockCGrQ" : "fPcg") + << " EXACT fp64" << std::endl; + } + + // One right-hand side through the LinearFunction interface (nrhs must be 1) + void Solve(const FineField &src, FineField &sol) + { + GRID_ASSERT( nrhs == 1 ); + sol.Checkerboard() = src.Checkerboard(); + GridStopWatch w; w.Start(); + Outer(src,sol); + w.Stop(); + FineField Ax(src.Grid()); Ax.Checkerboard()=src.Checkerboard(); + FineOp.HermOp(sol,Ax); + Ax=Ax-src; + std::cout << GridLogMessage << "HDCGSolver: single rhs total " << w.Elapsed() + << "; FINAL Nrhs 1 (LinearFunction): true residual = " << std::sqrt(norm2(Ax)/norm2(src)) << std::endl; + } + void Solve(std::vector &src, std::vector &sol) + { + GRID_ASSERT( (int)src.size() == nrhs ); + for(int r=0;r Shifted(HermOp,S.RefineShift); + HermOpAdaptor ShiftedOp(Shifted); + C.Coarsen(ShiftedOp,1); // refinement below still needs the RAW basis + C.CoarseLanczos(nr); + + HDCGMultiGridParams PR = C.Params; + PR.Outer.Tol = S.RefineTol; + PR.Outer.MaxIt = S.RefineMaxIt; + PR.Smoother.Order = S.RefineSmootherOrder; + HDCGSolver Refiner(Ddwf,DdwfF,C,PR,nr,S.RefineShift); + C.RefineSubspaceHDCG(HermOp,Refiner.Outer,nr); + + if ( S.RefinedSubspaceFile.length() ) + saveSubspace(C.rawNull, S.RefinedSubspaceFile, S.FileControl); +} + +NAMESPACE_END(Grid); diff --git a/Grid/algorithms/multigrid/HDCGMultiGridParams.h b/Grid/algorithms/multigrid/HDCGMultiGridParams.h new file mode 100644 index 000000000..2922050dc --- /dev/null +++ b/Grid/algorithms/multigrid/HDCGMultiGridParams.h @@ -0,0 +1,177 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./Grid/algorithms/multigrid/HDCGMultiGridParams.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); + +////////////////////////////////////////////////////////////////////////////////////// +// Parameters of the two-level mrhs HDCG (Hermitian, Schur-preconditioned) +// multigrid, HDCGMultiGrid.h. Serialisable: read from XML, template written +// on a missing file. Setup (blocking, batch, basis file, fine precision) is +// the same struct as the PVdagM chain's; Block2 is only used when a +// coarse-coarse level exists but must still divide the coarse lattice. +////////////////////////////////////////////////////////////////////////////////////// + +GRID_SERIALIZABLE_ENUM(HDCGSubspaceMethod, undef, ChebyshevNew, 1, Chebyshev, 2, FineEvecs, 3); +GRID_SERIALIZABLE_ENUM(HDCGRefineMethod, undef, None, 1, CG, 2, HDCG, 3); +GRID_SERIALIZABLE_ENUM(HDCGLanczosMethod, undef, None, 1, Block, 2, Plain, 3); +GRID_SERIALIZABLE_ENUM(HDCGCoarseMethod, undef, DeflatedCG, 1, Chebyshev, 2); +GRID_SERIALIZABLE_ENUM(HDCGSmootherType, undef, CG, 1, Chebyshev, 2); +GRID_SERIALIZABLE_ENUM(HDCGOuterAlgorithm, undef, fPcg, 1, PrecBlockCGrQ, 2); + +// Near-null basis. ChebyshevNew filters noise through the Hermitian operator +// with the upper bound only; Chebyshev takes lo/hi/order; FineEvecs sums +// FineEvecSample consecutive fp32 fine eigenvector records per basis vector. +// A ChebyHi <= 0 means 1.1 x the power-method estimate. FileControl is the +// scidac binary-IO control word of the basis files (0 for pre-2026 HDCG files). +struct HDCGSubspaceParams : Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(HDCGSubspaceParams, + HDCGSubspaceMethod, Method, + RealD, ChebyHi, + RealD, ChebyLo, + int, ChebyOrder, + std::string, FineEvecFile, + int, FineEvecSample, + int, FileControl, + HDCGRefineMethod, Refine, + RealD, RefineShift, // CG: the shift of the refining CG; HDCG: of the refining operator + RealD, RefineTol, + int, RefineMaxIt, + int, RefineSmootherOrder, // HDCG refine only + std::string, RefinedSubspaceFile); + HDCGSubspaceParams() + : Method(HDCGSubspaceMethod::ChebyshevNew), ChebyHi(0.0), ChebyLo(0.01), ChebyOrder(500), + FineEvecFile(""), FineEvecSample(1), + FileControl(BinaryIO::BINARYIO_LEXICOGRAPHIC|BinaryIO::BINARYIO_AGGREGATE), + Refine(HDCGRefineMethod::None), RefineShift(1.0e-3), RefineTol(1.0e-3), RefineMaxIt(500), + RefineSmootherOrder(11), RefinedSubspaceFile("") {}; +}; + +// Coarse eigenvectors for the deflated coarse solve: block Lanczos at the +// solve's Nrhs, or plain Lanczos at Nrhs 1, on the same coarse operator. +// ChebyHi <= 0 means 1.1 x the power-method estimate on the coarse operator. +struct HDCGLanczosParams : Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(HDCGLanczosParams, + HDCGLanczosMethod, Method, + int, Nstop, + int, Nk, + int, Nm, + RealD, ChebyLo, + RealD, ChebyHi, + int, ChebyOrder, + RealD, Tol, + int, MaxIt, + std::string, EvecFile); + HDCGLanczosParams() + : Method(HDCGLanczosMethod::Block), Nstop(32), Nk(32), Nm(128), + ChebyLo(0.01), ChebyHi(0.0), ChebyOrder(201), Tol(1.0e-4), MaxIt(100), EvecFile("") {}; +}; + +struct HDCGCoarseSolverParams : Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(HDCGCoarseSolverParams, + HDCGCoarseMethod, Method, + RealD, Tol, // DeflatedCG + int, MaxIt, + RealD, ChebyLo, // Chebyshev polynomial inverse + RealD, ChebyHi, + int, ChebyOrder); + HDCGCoarseSolverParams() + : Method(HDCGCoarseMethod::DeflatedCG), Tol(5.0e-2), MaxIt(3000), + ChebyLo(1.0e-2), ChebyHi(40.0), ChebyOrder(120) {}; +}; + +// CG: Order iterations of CG on the operator shifted by Shift. +// Chebyshev: order Order polynomial on [ChebyLo,ChebyHi] of the unshifted operator. +struct HDCGSmootherParams : Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(HDCGSmootherParams, + HDCGSmootherType, Type, + RealD, Shift, + int, Order, + RealD, ChebyLo, + RealD, ChebyHi); + HDCGSmootherParams() : Type(HDCGSmootherType::CG), Shift(2.0), Order(7), ChebyLo(2.0), ChebyHi(92.0) {}; +}; + +struct HDCGOuterParams : Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(HDCGOuterParams, + RealD, Tol, + int, MaxIt, + HDCGOuterAlgorithm, Algorithm); + HDCGOuterParams() : Tol(1.0e-8), MaxIt(500), Algorithm(HDCGOuterAlgorithm::fPcg) {}; +}; + +struct HDCGMultiGridParams : Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(HDCGMultiGridParams, + MGSetupParams, Setup, + HDCGSubspaceParams, Subspace, + HDCGLanczosParams, Lanczos, + HDCGCoarseSolverParams, CoarseSolver, + HDCGSmootherParams, Smoother, + HDCGOuterParams, Outer); + HDCGMultiGridParams() {}; +}; + +inline void CheckValidity(const HDCGMultiGridParams &P) +{ + GRID_ASSERT( P.Setup.Block1.size() == 4 ); + GRID_ASSERT( P.Setup.Block2.size() == 4 ); + GRID_ASSERT( P.Setup.CoarsenBatch >= 1 ); + GRID_ASSERT( P.Setup.FinePrecision != MGPrecision::undef ); + GRID_ASSERT( P.Subspace.Method != HDCGSubspaceMethod::undef ); + GRID_ASSERT( P.Subspace.Refine != HDCGRefineMethod::undef ); + GRID_ASSERT( P.Subspace.FineEvecSample >= 1 ); + GRID_ASSERT( P.Lanczos.Method != HDCGLanczosMethod::undef ); + GRID_ASSERT( P.Lanczos.Nstop <= P.Lanczos.Nk && P.Lanczos.Nk <= P.Lanczos.Nm ); + GRID_ASSERT( P.CoarseSolver.Method != HDCGCoarseMethod::undef ); + GRID_ASSERT( P.CoarseSolver.Tol > 0.0 ); + GRID_ASSERT( P.Smoother.Type != HDCGSmootherType::undef ); + GRID_ASSERT( P.Smoother.Order > 0 ); + GRID_ASSERT( P.Outer.Tol > 0.0 ); + GRID_ASSERT( P.Outer.Algorithm != HDCGOuterAlgorithm::undef ); +} + +// As ReadPVdagMMultiGridParams: a missing file gets a template written next +// to it and false returned; an empty filename means defaults. +inline bool ReadHDCGMultiGridParams(HDCGMultiGridParams &P, const std::string &xmlfile) +{ + if ( xmlfile.length() ) { + bool good; + { std::ifstream f(xmlfile); good = f.good(); } + if ( !good ) { + std::cout << GridLogMessage << "HDCGMultiGridParams: " << xmlfile << " does not exist" << std::endl; + if ( GlobalSharedMemory::WorldRank == 0 ) { + XmlWriter WR(xmlfile+".templ"); + write(WR, "HDCGMultiGridParams", P); + std::cout << GridLogMessage << "HDCGMultiGridParams: template written to " + << xmlfile << ".templ" << std::endl; + } + return false; + } + XmlReader RD(xmlfile); + read(RD, "HDCGMultiGridParams", P); + } + CheckValidity(P); + std::cout << GridLogMessage << "HDCGMultiGridParams (" + << (xmlfile.length() ? xmlfile : std::string("defaults")) << "):" << std::endl; + std::cout << P << std::endl; + return true; +} + +NAMESPACE_END(Grid); diff --git a/Grid/algorithms/multigrid/MrhsPreconditioner.h b/Grid/algorithms/multigrid/MrhsPreconditioner.h new file mode 100644 index 000000000..90121f5c6 --- /dev/null +++ b/Grid/algorithms/multigrid/MrhsPreconditioner.h @@ -0,0 +1,73 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./Grid/algorithms/multigrid/MrhsPreconditioner.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); + +////////////////////////////////////////////////////////////////////// +// The mrhs interfaces. This header is in the Algorithms.h include +// chain (AdefMrhs.h needs it); the concrete mrhs classes are in +// MrhsMultiGrid.h, reached through PVdagMMultiGrid.h / HDCGMultiGrid.h. +// +// MrhsLinearFunction: 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; + virtual ~MrhsLinearFunction(){}; +}; + +////////////////////////////////////////////////////////////////////// +// An mrhs preconditioner for the two-level CG family (AdefMrhs.h): +// operator() is M1; Vstart chooses x_0 from the source (Tang's ADEF-2 +// start x_0 = Q b; the default is the zero start). The timers are the +// M1 breakdown the outer solver reports. +////////////////////////////////////////////////////////////////////// +template +class MrhsPreconditioner : public MrhsLinearFunction { +public: + GridStopWatch ProjectTimer, PromoteTimer, DeflateTimer, CoarseTimer, FineTimer, SmoothTimer; + virtual void Vstart(std::vector &x, std::vector &src) + { + for(int r=0;r<(int)x.size();r++) x[r] = Zero(); + } + virtual void ResetTimers(void) + { + ProjectTimer.Reset(); PromoteTimer.Reset(); DeflateTimer.Reset(); + CoarseTimer.Reset(); FineTimer.Reset(); SmoothTimer.Reset(); + } + virtual void ReportTimers(const std::string &prefix) + { + std::cout< + + 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); + +////////////////////////////////////////////////////////////////////// +// Multigrid setup I/O shared by the PVdagM and HDCG chains: a set of +// bare-vector scidac records (a near-null basis, or a set of +// eigenvectors with their eigenvalues beside them in XML). +// +// The coarse OPERATOR is not written: coarsening is cheap relative +// to producing the basis, and is redone. +// +// `control` is the scidac binary-IO control word. The library default +// is lexicographic + aggregated; pass 0 to read files written by the +// pre-2026 HDCG drivers, which used control 0. +////////////////////////////////////////////////////////////////////// +static const int MultiGridIODefaultControl = BinaryIO::BINARYIO_LEXICOGRAPHIC|BinaryIO::BINARYIO_AGGREGATE; + +template +void saveSubspace(std::vector &subspace, std::string const fname, + int control = MultiGridIODefaultControl) +{ +#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, 0, control); + } + SW.close(); +#endif +} +template +void loadSubspace(std::vector &subspace, std::string const fname, + int control = MultiGridIODefaultControl) +{ +#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, control); + } + SR.close(); +#endif +} + +////////////////////////////////////////////////////////////////////// +// Eigenpairs: the vectors as a subspace file, the values in +// .evals.xml. evec.size() fixes how many are read. +////////////////////////////////////////////////////////////////////// +template +void saveEigenpairs(std::vector &evec, std::vector &eval, std::string const fname, + int control = MultiGridIODefaultControl) +{ + GRID_ASSERT(evec.size()==eval.size()); + saveSubspace(evec, fname, control); + if ( evec[0].Grid()->IsBoss() ) { + XmlWriter WR(fname+".evals.xml"); + write(WR,"evals",eval); + } +} +// The eigenvalue file sets the count; the vectors are allocated on grid. +template +void loadEigenpairs(std::vector &evec, std::vector &eval, std::string const fname, + GridBase *grid, int control = MultiGridIODefaultControl) +{ + XmlReader RD(fname+".evals.xml"); + read(RD,"evals",eval); + evec.clear(); + evec.resize(eval.size(), Field(grid)); + loadSubspace(evec, fname, control); +} + +NAMESPACE_END(Grid);