From 97c9178785571cd46e4e495e00e0dc99f04dafc7 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Thu, 20 Aug 2026 11:44:26 -0400 Subject: [PATCH] Ready for faster 3 level solve with dense coarse and BLAS on SRHS ! Exciting --- CLAUDE.md | 2 + Grid/algorithms/multigrid/DenseCoarseMatrix.h | 160 +++++-- .../multigrid/GeneralCoarsenedMatrix.h | 11 + .../GeneralCoarsenedMatrixMultiRHS.h | 15 +- .../GeneralCoarsenedMatrixMultiRHSV2.h | 15 +- ...le_pvdagm_mrhs_3level_DenseCoarseMatrix.cc | 5 +- ...mple_pvdagm_v2_3level_DenseCoarseMatrix.cc | 411 +++++++++++++++++- tests/debug/Test_coarse_v2_coarsen.cc | 31 ++ 8 files changed, 573 insertions(+), 77 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index 941667a86..741262340 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -174,6 +174,8 @@ Every programme is wrapped in `Grid_init(&argc, &argv)` / `Grid_finalize()` (`Gr - Template structure: most classes are templated on `<_FImpl>` (fermion impl) or `` (gauge impl), which encode the representation and precision. Instantiation is controlled by `--enable-fermion-instantiations`. - **Tensor indices are positional, not labelled.** The `Grid/tensors/` arithmetic recurses structurally over the `iScalar`/`iVector`/`iMatrix` nest: each level defines only the {scalar,vector,matrix}² products at its own level, with element types resolved by automatic type deduction, so every colour/spin/lorentz combination composes from ~200 lines (versus the pre-C++11 QDP++/PETE approach of machine-generating every case). An index's meaning derives entirely from its nesting depth counted from the outside; `iScalar` is the identity/broadcast case at every level. Never insert or remove a nesting level casually — the multiplication tables contract by position. - **Multigrid coarsening deepens the tensor nest by one level.** A coarse site vector is `iVector`, and `innerProduct` on it returns `iScalar` — one level deeper than the fine block scalar. So the block-inner-product scalar type gains one `iScalar` wrapper per MG level (fine: `vTComplex`; level 2: `iScalar`; see `examples/Example_pvdagm_3level.cc`). When calling `blockInnerProduct`/`blockZAXPY`/`blockOrthogonalise` on coarse fields, the coarse scalar type must match `decltype(innerProduct(siteVector(),siteVector()))` exactly; a wrong depth fails to compile (no viable `operator=` deep in the instantiation chain) rather than mis-contracting. +- **Grids are borrowed, never owned.** `conformable` is pointer identity, so every object that interoperates must hold the *same* `GridCartesian *`; a class that minted its own grid internally could never conform with anything else. Ownership is therefore not available, and lifetime is managed by scope discipline instead of reference counting: whoever creates a grid retains it beyond every object it handed a reference to. Anything *derived* from a grid inherits this — `~PaddedCell` dereferences its `unpadded_grid`, so a `PaddedCell` cannot even be **destroyed** after its parent grid, only used. Where a consumer must let go early, it offers an explicit hand-back (`MultiGeneralCoarsenedOperatorV2::ReleaseGrid()`) to be called *before* the grid is destroyed. + - The `RealD`/`RealF`/`ComplexD`/`ComplexF` typedefs are used everywhere; avoid raw `double`/`float`. - Use `GRID_ASSERT(cond)` (defined in `Grid/GridStd.h`), not bare `assert` — it prints a Grid-formatted message and aborts cleanly under MPI. - Logging is stream-based, not macro-based: `std::cout << GridLogMessage << ... << std::endl;`. Channels declared in `Grid/log/Log.h` include `GridLogError`, `GridLogWarning`, `GridLogDebug`, `GridLogPerformance`, `GridLogIterative`, `GridLogSolver`, `GridLogHMC`, `GridLogComms`, `GridLogMemory`, `GridLogDslash`, `GridLogIRL`, `GridLogMG`. A subset is switched on at runtime with e.g. `--log Error,Warning,Message,Performance,Iterative,Integrator,Debug,Colours` (names given without the `GridLog` prefix). diff --git a/Grid/algorithms/multigrid/DenseCoarseMatrix.h b/Grid/algorithms/multigrid/DenseCoarseMatrix.h index 174a270e8..780f508ca 100644 --- a/Grid/algorithms/multigrid/DenseCoarseMatrix.h +++ b/Grid/algorithms/multigrid/DenseCoarseMatrix.h @@ -81,19 +81,25 @@ NAMESPACE_BEGIN(Grid); // Tensor-depth agnostic: site scalar objects treated as contiguous ComplexD // (iScalar wrappers add no data), so any MG level's coarse operator imports. ////////////////////////////////////////////////////////////////////////////////////// -template -class DenseCoarseMatrix : public LinearFunction::CoarseVector> { +// +// Depends on a coarse operator only to extract its matrix elements; thereafter +// it is given a coarse vector and applies the inverse. Import() is a template +// member so any of the coarse classes will do, and the type of a +// DenseCoarseMatrix does not record which one built it. +// +template +class DenseCoarseMatrix : public LinearFunction > > { public: - typedef GeneralCoarsenedMatrix GeneralCoarseOp; - typedef typename GeneralCoarseOp::CoarseVector Field; - typedef typename GeneralCoarseOp::CoarseMatrix CoarseMatrix; + typedef iVector siteVector; + typedef Lattice CoarseVector; + typedef Lattice > CoarseMatrix; + typedef CoarseVector Field; using LinearFunction::operator(); typedef typename Field::vector_object vobj; typedef typename vobj::scalar_object sobj; typedef typename CoarseMatrix::vector_object Mvobj; typedef typename Mvobj::scalar_object Msobj; - GeneralCoarseOp &_Op; // the coarse operator: stencil source + certificate oracle GridBase *grid; int nd; int64_t N; // dense rank = gSites * nbasis @@ -121,12 +127,11 @@ public: int devSum; double schurAuditRel; // DENSE_SCHUR=2: rel slab diff single-vs-schur (-1 = not run) - DenseCoarseMatrix(GeneralCoarseOp &Op, GridBase *g) - : _Op(Op), grid(g) + DenseCoarseMatrix(GridBase *g) + : grid(g) { GRID_ASSERT( sizeof(sobj) == nbasis*sizeof(ComplexD) ); GRID_ASSERT( sizeof(Msobj) == nbasis*nbasis*sizeof(ComplexD) ); - GRID_ASSERT( grid == Op.Grid() ); nd = grid->_ndimension; N = grid->gSites() * nbasis; lsites = grid->lSites(); @@ -158,7 +163,16 @@ public: } slab.resize((uint64_t)nrows * N); + } + //////////////////////////////////////////////////////////////////// + // The only place a coarse operator is needed: pull its elements, invert, + // and make the slab resident. Any class exposing Geometry() and + // ExtractMatrix(p,A) will do -- single RHS or either multiRHS. + //////////////////////////////////////////////////////////////////// + template + void Import(CoarseOp &Op) + { double t0 = usecond(); //////////////////////////////////////////////////////////////////// // 0. Slab cache: SLAB_FILE= -> per-rank raw file .. @@ -187,9 +201,9 @@ public: } } if (!loaded) { - ImportDense(); // slab <- my rows of A (LOCAL, no comms) - ImportCertificate(); // dense apply == Op.M on a non-constant vector - InvertDense(); // slab <- my rows of A^{-1} + ImportDense(Op); // slab <- my rows of A (LOCAL, no comms) + ImportCertificate(Op); // dense apply == Op.M, before inversion + InvertDense(Op); // slab <- my rows of A^{-1} double t1 = usecond(); if (sfile) { FILE *f = fopen(slabfile.c_str(),"wb"); @@ -251,7 +265,7 @@ public: double ta = usecond(); (*this)(x, y); double tb = usecond(); - _Op.M(y, z); + ApplyOracle(Op, y, z); z = z - x; RealD rel = std::sqrt(norm2(z)/norm2(x)); std::cout << GridLogMessage << "DenseCoarseMatrix: VERIFY ||A Ainv x - x||/||x|| = " @@ -262,11 +276,52 @@ public: << (usecond()-t0)/1.0e6 << " s" << std::endl; } + //////////////////////////////////////////////////////////////////// + // Apply the source operator to a D dimensional field, whichever kind it is. + // + // A multiRHS operator lives on the D+1 grid, so drive it with several right + // hand sides at once: slice r carries (r+1)*in, and linearity says the + // results must scale likewise. One apply, and unlike a single rhs check it + // also catches rhs mixing. Cheap check, not a production path. + //////////////////////////////////////////////////////////////////// + template + void ApplyOracle(CoarseOp &Op,const Field &in, Field &out) + { + if ( Op.Grid() == grid ) { Op.M(in,out); return; } + + GridBase *mgrid = Op.Grid(); + GRID_ASSERT(mgrid->_ndimension == nd+1); + int nr = mgrid->_fdimensions[0]; + + Field min(mgrid), mout(mgrid); + for(int r=0;r= 1.0e-6 ) { + std::cout << GridLogMessage << "DenseCoarseMatrix: oracle rhs "< dense import of MY ROWS of A (no comms): // Dense[(s,a),(wrap(s+shift_p),b)] += A[p][s]_{a,b} //////////////////////////////////////////////////////////////////// - void ImportDense(void) + template + void ImportDense(CoarseOp &Op) { double t = -usecond(); Coordinate gdims = grid->GlobalDimensions(); @@ -276,12 +331,12 @@ public: uint64_t nelem = (uint64_t)nrows * N; thread_for(i, nelem, { slab[i] = ComplexF(0.0,0.0); }); - for(int p=0; p<_Op.geom.npoint; p++){ - Coordinate shift = _Op.geom.shifts[p]; + for(int p=0; pGlobalSumVector(&gpk, 1); std::cout << GridLogMessage << "DenseCoarseMatrix: DEBUG p=" << p - << " norm2(_A[p]) " << norm2(_Op._A[p]) + << " norm2(_A[p]) " << norm2(Aun) << " norm2(Extract) " << norm2(Aun) << " sum|peek|^2 " << gpk << std::endl; } @@ -346,7 +401,7 @@ public: grid->GlobalSumVector(&gz, 1); std::cout << GridLogMessage << "DenseCoarseMatrix: stencil->dense import took " - << t/1.0e6 << " s (" << _Op.geom.npoint << " points, local, no comms)" + << t/1.0e6 << " s (" << Op.Geometry().npoint << " points, local, no comms)" << " zero rows " << (int64_t)gz << "/" << N << std::endl; // Debug: coordinate pattern of live sites (mechanism fingerprint) @@ -375,7 +430,8 @@ public: // 2. IMPORT CERTIFICATE: dense rows vs Op.M on a NON-CONSTANT vector. // (Constant x has x[s+d]==x[s-d]: blind to a shift-sign error.) //////////////////////////////////////////////////////////////////// - void ImportCertificate(void) + template + void ImportCertificate(CoarseOp &Op) { Field x(grid); Field Ax(grid); Field Dx(grid); for(int ss=0; ss ref(slab); // Ainv, single path slab = Aimp; - InvertDenseSchur(); // slab = Ainv, Schur path + InvertDenseSchur(Op); // slab = Ainv, Schur path // NaN-PROOF comparison: max() masks NaN, so count non-finite // entries in each result explicitly. @@ -713,17 +770,18 @@ public: // different precision order). NaN-proof: non-finite entries are // counted explicitly since max() silently masks NaN. //////////////////////////////////////////////////////////////////// - void ImportDenseFP64(BlockRows &S, std::vector &g2rm) + template + void ImportDenseFP64(CoarseOp &Op, BlockRows &S, std::vector &g2rm) { Coordinate gdims = grid->GlobalDimensions(); int sign = getenv("DENSE_IMPORT_SIGN") ? atoi(getenv("DENSE_IMPORT_SIGN")) : 1; GRID_ASSERT( sign==1 || sign==-1 ); std::vector h((uint64_t)nrows*N, ComplexD(0.0,0.0)); - for(int p=0; p<_Op.geom.npoint; p++) + for(int p=0; p + void InvertDenseSchur(CoarseOp &Op) { double t1 = usecond(); int P = grid->ProcessorCount(); @@ -812,7 +871,7 @@ public: } BlockRows S; - ImportDenseFP64(S, g2rm); + ImportDenseFP64(Op, S, g2rm); int64_t panelBytes = getenv("DENSE_PANEL_BYTES") ? atol(getenv("DENSE_PANEL_BYTES")) : (int64_t)1024*1024*1024; @@ -921,13 +980,21 @@ public: ((ComplexD *)&s)[b] = ComplexD(hY[ss*nbasis + b]); pokeLocalSite(s, psi, myLcoor[ss]); } - if ( getenv("DENSE_CC_CHECK") ) { - Field tmp(grid); - _Op.M(psi, tmp); - tmp = tmp - src; - std::cout << GridLogMessage << "DenseCoarseMatrix: apply defect ||A x - b||/||b|| = " - << std::sqrt(norm2(tmp)/norm2(src)) << std::endl; - } + } + + //////////////////////////////////////////////////////////////////// + // Defect of an applied inverse. Was a DENSE_CC_CHECK block inside + // operator(), but that is virtual and cannot take an operator, so the + // caller now asks for it explicitly. + //////////////////////////////////////////////////////////////////// + template + void CheckApply(CoarseOp &Op,const Field &src,const Field &psi) + { + Field tmp(grid); + Op.M(psi, tmp); + tmp = tmp - src; + std::cout << GridLogMessage << "DenseCoarseMatrix: apply defect ||A x - b||/||b|| = " + << std::sqrt(norm2(tmp)/norm2(src)) << std::endl; } //////////////////////////////////////////////////////////////////// @@ -960,14 +1027,17 @@ public: double t1 = usecond(); std::cout << GridLogMessage << "DenseCoarseMatrix: batched apply " << nr << " rhs took " << (t1-t0)/1000.0 << " ms (" << (t1-t0)/1000.0/nr << " ms/rhs)" << std::endl; - if ( getenv("DENSE_CC_CHECK") ) { - Field tmp(grid); - for(int rr=0; rrImport(LittleDiracOpL2); MrhsDenseCC.reset(new MrhsDenseCCSolve(*DenseCC, CoarseCoarse5d, nrhs)); } diff --git a/examples/Example_pvdagm_v2_3level_DenseCoarseMatrix.cc b/examples/Example_pvdagm_v2_3level_DenseCoarseMatrix.cc index 6df7929f1..01ea463cb 100644 --- a/examples/Example_pvdagm_v2_3level_DenseCoarseMatrix.cc +++ b/examples/Example_pvdagm_v2_3level_DenseCoarseMatrix.cc @@ -29,8 +29,8 @@ Author: Peter Boyle // // PVdagM three level multigrid on the V2 coarse operator. // -// STAGE ONE: grids, types, subspace, and the L1 coarsening only. The L2 -// chain, the dense bottom and the solves are not here yet. +// STAGES ONE AND TWO: grids, types, subspace, and the L1 and L2 coarsenings. +// The dense bottom and the solves are not here yet. // // Differences from Example_pvdagm_mrhs_3level_DenseCoarseMatrix.cc: // @@ -53,14 +53,15 @@ Author: Peter Boyle // aggregation gives e_k, and the near null content is silently gone. The // || - I||_F guard below is what catches that. // -// Env: LATT LS MASS NBASIS(compile time) NRHS BLOCK COARSEN_BATCH -// HOT_START CONFIG SUBSPACE_FILE V1_CHECK +// Env: LATT LS MASS NBASIS(compile time) NRHS BLOCK BLOCK2 COARSEN_BATCH +// HOT_START CONFIG SUBSPACE_FILE V1_CHECK MRHS_COARSEN // #include #include #include #include +#include #include @@ -78,12 +79,32 @@ int Ls = 24; int CoarsenBatch = 9; std::vector lat_size({48,48,48,96}); +// Solver tuning, values as in the V1 example +RealD FineSmootherShift = 0.1; +int FineSmootherOrder = 16; +RealD CoarseSmootherShift = 0.1; +int CoarseSmootherNstep = 4; +RealD CoarseSolverTol = 0.03; +int CoarseSolverOrder = 200; +RealD OuterTol = 1.0e-8; +int OuterMmax = 8; +int OuterNstep = 8; + void ParseEnvironment(void) { if(getenv("MASS")) mass = atof(getenv("MASS")); if(getenv("NRHS")) Nrhs = atoi(getenv("NRHS")); if(getenv("LS")) Ls = atoi(getenv("LS")); if(getenv("COARSEN_BATCH")) CoarsenBatch= atoi(getenv("COARSEN_BATCH")); + if(getenv("FineSmootherShift")) FineSmootherShift = atof(getenv("FineSmootherShift")); + if(getenv("FineSmootherOrder")) FineSmootherOrder = atoi(getenv("FineSmootherOrder")); + if(getenv("CoarseSmootherShift"))CoarseSmootherShift= atof(getenv("CoarseSmootherShift")); + if(getenv("CoarseSmootherNstep"))CoarseSmootherNstep= atoi(getenv("CoarseSmootherNstep")); + if(getenv("CoarseSolverTol")) CoarseSolverTol = atof(getenv("CoarseSolverTol")); + if(getenv("CoarseSolverOrder")) CoarseSolverOrder = atoi(getenv("CoarseSolverOrder")); + if(getenv("OuterTol")) OuterTol = atof(getenv("OuterTol")); + if(getenv("OuterMmax")) OuterMmax = atoi(getenv("OuterMmax")); + if(getenv("OuterNstep")) OuterNstep = atoi(getenv("OuterNstep")); if(getenv("LATT")){ Coordinate l; GridCmdOptionIntVector(std::string(getenv("LATT")),l); @@ -139,9 +160,8 @@ public: }; ////////////////////////////////////////////////////////////////////// -// || - I||_F over a set of coarse vectors. ~0.23 means the raw near -// null content survived the projection; ~sqrt(N_sites) means block -// orthonormal vectors leaked in and every image collapsed to e_k. +// || - I||_F over a set of coarse vectors. Small means the raw near null +// content survived the projection; see GramGuard for where a leak lands. ////////////////////////////////////////////////////////////////////// template RealD GramDefect(std::vector &v) @@ -157,6 +177,228 @@ RealD GramDefect(std::vector &v) 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 ); +} + + +////////////////////////////////////////////////////////////////////// +// Shifted variants for the smoothers +////////////////////////////////////////////////////////////////////// +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) { 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 Op (const Field &in, Field &out){ Field tmp(in.Grid()); _Mat.M(in,tmp); _PV.Mdag(tmp,out); out = out + shift*in; } + void AdjOp (const Field &in, Field &out){ Field tmp(in.Grid()); _PV.M(tmp,out); _Mat.Mdag(in,tmp); out = out + shift*in; } + void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ assert(0); } + void HermOp(const Field &in, Field &out){ Field tmp(in.Grid()); Op(in,tmp); AdjOp(tmp,out); } +}; + +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) { 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 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){ assert(0); } + void HermOp (const Field &in, Field &out) { Field tmp(in.Grid()); Op(in,tmp); AdjOp(tmp,out); } +}; + +////////////////////////////////////////////////////////////////////// +// Dense L3 solve on the packed D+1 coarse-coarse field +////////////////////////////////////////////////////////////////////// +template +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); + } +}; + +////////////////////////////////////////////////////////////////////// +// mrhs interfaces + single-polynomial mrhs PGCR +////////////////////////////////////////////////////////////////////// +template +class MrhsLinearFunction { +public: + virtual void operator()(std::vector &in, std::vector &out) = 0; +}; + +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; + 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){ 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,b,rq; RealD zAAz; int nrhs=src.size(); GridBase *grid=src[0].Grid(); + std::vector r(nrhs,grid),z(nrhs,grid),Az(nrhs,grid); + 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); + for(int back=0;back=0); + b=-real(vinnerProduct(q[peri_back],Az))/qq[peri_back]; + vaxpy(p[peri_kp],b,p[peri_back],p[peri_kp]); vaxpy(q[peri_kp],b,q[peri_back],q[peri_kp]); } + qq[peri_kp]=vnorm2(q[peri_kp]); + } + GRID_ASSERT(0); return cp; + } +}; + +////////////////////////////////////////////////////////////////////// +// L2->L3 mrhs V-cycle on the D+1 coarse field +////////////////////////////////////////////////////////////////////// +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) { + int nrhs=_nrhs; + 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); + + CCsol=Zero(); + _CoarseCoarseSolve(CCsrc,CCsol); + + _Projector.blockPromote(vec1,CCsol); + add(out,out,vec1); + + _CoarseOp.Op(out,vec1); sub(vec1,in,vec1); + vec2=Zero(); + _CoarseSmoother(vec1,vec2); + add(out,out,vec2); + } +}; + +////////////////////////////////////////////////////////////////////// +// L1->L2 mrhs V-cycle +////////////////////////////////////////////////////////////////////// +template +class MrhsTwoLevelMG : public MrhsLinearFunction { +public: + typedef MrhsCoarseVector CoarseVector; + LinearOperatorBase &_FineOperator; + FineSmoother &_PostSmoother; + MultiRHSBlockProject &_Projector; + LinearFunction &_CoarseSolve; + GridBase *_CoarseGrid, *_CoarseGridMrhs; + 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){ + 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); + _Projector.blockProject(vec1,CsrcMrhs); + + CsolMrhs=Zero(); + _CoarseSolve(CsrcMrhs,CsolMrhs); + + _Projector.blockPromote(vec1,CsolMrhs); + for(int r=0;r PVdagM_t; + typedef PVdagMLinearOperator PVdagM_t; + typedef ShiftedPVdagMLinearOperator ShiftedPVdagM_t; PVdagM_t PVdagM(Ddwf,Dpv); ////////////////////////////////////////////////////////////////////// @@ -310,13 +553,7 @@ int main (int argc, char ** argv) MrhsProjector.blockProject(rawNull,psi_coarse); // RAW vectors in rawNull.clear(); rawNull.shrink_to_fit(); - { - RealD defect = GramDefect(psi_coarse); - RealD leak = std::sqrt((double)Coarse5d->gSites()); - std::cout << GridLogMessage << "GUARD: || - I||_F = " << defect - << " (~0.23 good; ~sqrt(N_coarse)=" << leak << " = e_k leak)" << std::endl; - GRID_ASSERT( defect < leak ); - } + GramGuard("psi_coarse",psi_coarse,Coarse5d); ////////////////////////////////////////////////////////////////////// // Optional cross check of the coarse matrix elements against the V1 @@ -366,9 +603,9 @@ int main (int argc, char ** argv) ComplexD *w2=(ComplexD *)&h2[0]; int64_t words = sites*sizeof(calcMatrix)/sizeof(ComplexD); for(int64_t i=0;i - I||_F = " << defect - << " (~0.23 good; ~sqrt(N_cc)=" << leak << " = e_k leak)" << std::endl; - GRID_ASSERT( defect < leak ); + GramGuard("psi_cc",psi_cc,CoarseCoarse5d); } rawPsi.clear(); rawPsi.shrink_to_fit(); @@ -470,7 +703,137 @@ int main (int argc, char ** argv) GRID_ASSERT( norm2(ccout) > 0.0 ); } - std::cout << GridLogMessage << "*** stage two complete: L1 and L2 coarse operators built ***" << std::endl; + ////////////////////////////////////////////////////////////////////// + // STAGE THREE (part one): the dense bottom on L2. + // + // DenseCoarseMatrix is bilingual: it takes the elements through + // Geometry()/ExtractMatrix(), so the V2 operator serves directly. It does + // detect that a multiRHS op cannot apply on the D dimensional grid and + // skips its own certificate and VERIFY, so the equivalent check is done + // here instead, driving the L2 operator at Nrhs 1 through a slice. + ////////////////////////////////////////////////////////////////////// + typedef DenseCoarseMatrix DenseCC_t; + std::unique_ptr DenseCC; + + if ( getenv("DENSE_CC")==nullptr || atoi(getenv("DENSE_CC")) ) { + + std::cout << GridLogMessage << "*** L3 dense bottom: import from the V2 L2 operator ***" << std::endl; + DenseCC.reset(new DenseCC_t(CoarseCoarse5d)); + DenseCC->Import(CoarseOpL2); + + //////////////////////////////////////////////////////////////////// + // ||A Ainv x - x|| / ||x||, the check Import could not run itself + //////////////////////////////////////////////////////////////////// + Coordinate cc1latt({1,1,cclatt[0],cclatt[1],cclatt[2],cclatt[3]}); + GridCartesian *CoarseCoarseOne = new GridCartesian(cc1latt,cmsimd,cmmpi); + + CoarseOpL2.SetGrid(CoarseCoarseOne); + + CoarseCoarseVector x(CoarseCoarse5d),y(CoarseCoarse5d),z(CoarseCoarse5d); + GridParallelRNG dRNG(CoarseCoarse5d); dRNG.SeedFixedIntegers({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 << "L3 dense: ||A Ainv x - x||/||x|| = " << rel << std::endl; + GRID_ASSERT( rel < 1.0e-2 ); + + CoarseOpL2.SetGrid(CoarseCoarseMrhs); + delete CoarseCoarseOne; + } + + ////////////////////////////////////////////////////////////////////// + // STAGE THREE (part two): the solves. + // + // Both operators are driven from the SAME objects at whatever Nrhs is + // asked for -- the matrix elements were built once and survive SetGrid -- + // so single RHS and multiRHS are the same code path with a different grid. + ////////////////////////////////////////////////////////////////////// + GRID_ASSERT(DenseCC != nullptr); // the PGCR bottom is not ported yet + + typedef PrecGeneralisedConjugateResidualNonHermitian FineSmoother_t; + + ShiftedPVdagM_t ShiftedPVdagM(FineSmootherShift,Ddwf,Dpv); + TrivialPrecon simple_fine; + TrivialPrecon simpleC; + + auto RunSolve = [&](int nr) + { + std::cout << GridLogMessage << "**********************************************" << std::endl; + std::cout << GridLogMessage << " V2 THREE-level solve, Nrhs = " << nr << std::endl; + std::cout << GridLogMessage << "**********************************************" << std::endl; + + Coordinate cml({nr,1,clatt[0],clatt[1],clatt[2],clatt[3]}); + Coordinate ccml({nr,1,cclatt[0],cclatt[1],cclatt[2],cclatt[3]}); + GridCartesian *CMrhs = new GridCartesian(cml, cmsimd,cmmpi); + GridCartesian *CCMrhs = new GridCartesian(ccml,cmsimd,cmmpi); + + CoarseOpPV.SetGrid(CMrhs); + CoarseOpL2.SetGrid(CCMrhs); + + NonHermitianLinearOperator LinOpC (CoarseOpPV); + NonHermitianLinearOperator LinOpCC(CoarseOpL2); + + MrhsDenseCCSolve ccSolve(*DenseCC,nr); + + ShiftedLinearOperator ShiftedC(CoarseSmootherShift, LinOpC); + PrecGeneralisedConjugateResidualNonHermitian + CoarseSmootherGCR(0.01,1,ShiftedC,simpleC,CoarseSmootherNstep,CoarseSmootherNstep); + CoarseSmootherGCR.Level(2); CoarseSmootherGCR.Name("Csmoother"); CoarseSmootherGCR.SetZeroGuess(1); + + MrhsCoarseThreeLevelPrec + L2to3Precon(LinOpC, CoarseSmootherGCR, MrhsProjectorL2, ccSolve, + Coarse5d, CoarseCoarse5d, CCMrhs, nr); + + PrecGeneralisedConjugateResidualNonHermitian + L2PGCR(CoarseSolverTol, CoarseSolverOrder/16, LinOpC, L2to3Precon, 16, 16); + L2PGCR.Level(2); L2PGCR.Name("Couter"); L2PGCR.SetZeroGuess(1); + + FineSmoother_t SmootherGCR(0.0,1,ShiftedPVdagM,simple_fine,FineSmootherOrder,FineSmootherOrder); + SmootherGCR.Level(1); SmootherGCR.Name("Fsmoother"); SmootherGCR.SetZeroGuess(1); + + MrhsTwoLevelMG + ThreeLevelPrecon(PVdagM, SmootherGCR, MrhsProjector, L2PGCR, Coarse5d, CMrhs); + + MrhsPGCRNonHermitian + L1PGCR(OuterTol,1000,PVdagM,ThreeLevelPrecon,OuterMmax,OuterNstep); + L1PGCR.Level(1); L1PGCR.Name("Fouter"); L1PGCR.SetZeroGuess(1); + + std::vector src(nr,FGrid), sol(nr,FGrid); + for(int r=0;r