From 2900ce33b5630c61399df5796c2030e948f00695 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Wed, 19 Aug 2026 20:05:48 -0400 Subject: [PATCH] New files, including v2 multiRHS coarse op --- examples/Example_mdagm.cc | 240 +++++++ examples/Example_pvdagm_mrhs.cc | 662 ++++++++++++++++++ ...mple_pvdagm_v2_3level_DenseCoarseMatrix.cc | 476 +++++++++++++ 3 files changed, 1378 insertions(+) create mode 100644 examples/Example_mdagm.cc create mode 100644 examples/Example_pvdagm_mrhs.cc create mode 100644 examples/Example_pvdagm_v2_3level_DenseCoarseMatrix.cc diff --git a/examples/Example_mdagm.cc b/examples/Example_mdagm.cc new file mode 100644 index 000000000..aff55b316 --- /dev/null +++ b/examples/Example_mdagm.cc @@ -0,0 +1,240 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./examples/Example_mdagm.cc + + Copyright (C) 2023 + +Author: Peter Boyle + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution directory + *************************************************************************************/ + /* END LEGAL */ +#include +#include +#include + +using namespace std; +using namespace Grid; + +// Routes Op/AdjOp -> HermOp so that CoarsenOperator and CreateSubspace +// both see the HPD operator M†M rather than bare M. +template +class HermOpAdaptor : public LinearOperatorBase +{ + LinearOperatorBase &wrapped; +public: + HermOpAdaptor(LinearOperatorBase &wrapme) : wrapped(wrapme) {}; + void Op (const Field &in, Field &out) { wrapped.HermOp(in,out); } + void HermOp (const Field &in, Field &out) { wrapped.HermOp(in,out); } + void AdjOp (const Field &in, Field &out) { wrapped.HermOp(in,out); } + void OpDiag (const Field &in, Field &out) { 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 HermOpAndNorm(const Field &in, Field &out, RealD &n1, RealD &n2) { + wrapped.HermOp(in, out); + ComplexD dot = innerProduct(in, out); + n1 = real(dot); + n2 = norm2(out); + } +}; + +// Fixed-iteration CG smoother: runs exactly `iters` steps of CG on the +// shifted operator. tolerance=0 so CG never exits early. +template +class CGSmoother : public LinearFunction +{ +public: + using LinearFunction::operator(); + typedef LinearOperatorBase FineOperator; + FineOperator &_SmootherOperator; + int iters; + CGSmoother(int _iters, FineOperator &SmootherOperator) + : _SmootherOperator(SmootherOperator), iters(_iters) + { + std::cout << GridLogMessage << " CGSmoother order " << iters << std::endl; + } + void operator()(const Field &in, Field &out) + { + ConjugateGradient CG(0.0, iters, false); + out = Zero(); + CG(_SmootherOperator, in, out); + } +}; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + const int Ls = 24; + const int nbasis = 60; + const int cb = 0; + RealD M5 = 1.8; + RealD b = 1.5; + RealD c = 0.5; + + RealD mass = 0.00078; + { const char *e = getenv("MASS"); if (e && *e) mass = atof(e); } + + std::cout << GridLogMessage << "Mass: " << mass + << " Ls: " << Ls << " b=" << b << " c=" << c << std::endl; + std::cout << GridLogMessage << "nbasis: " << nbasis << std::endl; + + // ── Grids ────────────────────────────────────────────────────────────── + std::vector lat_size{48,48,48,96}; + + GridCartesian *UGrid = SpaceTimeGrid::makeFourDimGrid(lat_size, + GridDefaultSimd(Nd,vComplex::Nsimd()), + GridDefaultMpi()); + GridRedBlackCartesian *UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + GridCartesian *FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); + GridRedBlackCartesian *FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); + + Coordinate Block({4,4,3,4}); + Coordinate clatt = lat_size; + for (int d = 0; d < (int)clatt.size(); d++) clatt[d] /= Block[d]; + + GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, + GridDefaultSimd(Nd,vComplex::Nsimd()), + GridDefaultMpi()); + GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d); + + // ── RNGs ─────────────────────────────────────────────────────────────── + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers({5,6,7,8}); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers({1,2,3,4}); + + // ── Gauge field ──────────────────────────────────────────────────────── + LatticeGaugeField Umu(UGrid); + FieldMetaData header; + NerscIO::readConfiguration(Umu, header, std::string("/ccs/home/poare/ckpoint_lat.1000")); + + // ── Fermion operator ─────────────────────────────────────────────────── + MobiusFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c); + + MdagMLinearOperator MdagMOp(Ddwf); + HermOpAdaptor HermFineOp(MdagMOp); + + // ── Coarse geometry ──────────────────────────────────────────────────── + typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef LittleDiracOperator::CoarseVector CoarseVector; + typedef Aggregation Subspace; + + NextToNearestStencilGeometry5D geom(Coarse5d); + + // ── Power method: estimate upper end of M†M spectrum ─────────────────── + LatticeFermionD pm_src(FGrid); random(RNG5, pm_src); + PowerMethod PM; + RealD hi = PM(HermFineOp, pm_src); + std::cout << GridLogMessage << "Power method: hi = " << hi << std::endl; + + // ── Smoother: fixed-iteration CG on (M†M + lo) ───────────────────────── + // lo/hi ~ 2/95 matches the HDCG ratio; tune empirically. + RealD lo = hi / 40.0; + int ord = 12; + std::cout << GridLogMessage << "Smoother shift lo = " << lo + << " order = " << ord << std::endl; + ShiftedHermOpLinearOperator ShiftedFineOp(HermFineOp, lo); + CGSmoother Smoother(ord, ShiftedFineOp); + + // ── Subspace via CG inverse iteration ────────────────────────────────── + Subspace Aggregates(Coarse5d, FGrid, cb); + Aggregates.CreateSubspace(RNG5, HermFineOp, nbasis); + + // ── Cheap coarse deflation: diagonalise W BEFORE block-GS ────────────── + // Orthogonalise() applies block-local Gram-Schmidt which rotates subspace[i] + // into orthonormal block-local combinations, destroying the near-null + // property of individual vectors. We must extract the coarse zero-mode + // combinations χₖ = Σᵢ V[i,k] ψᵢ from the pre-GS near-null vectors first. + std::vector chi(nbasis, FGrid); + { + std::vector &psi = Aggregates.subspace; + LatticeFermionD tmp(FGrid); + Eigen::MatrixXcd W = Eigen::MatrixXcd::Zero(nbasis, nbasis); + for (int i = 0; i < nbasis; i++) { + HermFineOp.Op(psi[i], tmp); + for (int j = 0; j < nbasis; j++) + W(j, i) = TensorRemove(innerProduct(psi[j], tmp)); + } + Eigen::SelfAdjointEigenSolver esolver(W); + for (int k = 0; k < nbasis; k++) { + chi[k] = Zero(); + for (int i = 0; i < nbasis; i++) + chi[k] += ComplexD(esolver.eigenvectors()(i, k)) * psi[i]; + } + } + + Aggregates.Orthogonalise(); + + // ── Coarse operator ──────────────────────────────────────────────────── + LittleDiracOperator LittleDiracOp(geom, FGrid, Coarse5d); + LittleDiracOp.CoarsenOperator(HermFineOp, Aggregates); + + // ── Coarse linear operator ───────────────────────────────────────────── + HermitianLinearOperator LinOpCoarse(LittleDiracOp); + + // ── Project χₖ to coarse grid; Rayleigh quotients give deflation evals ─ + // ProjectToSubspace uses the post-GS basis, correctly mapping the pre-GS + // near-null combinations to coarse vectors via the block-local U†(x_c). + std::vector coarse_deflation_vecs(nbasis, Coarse5d); + std::vector coarse_deflation_evals(nbasis); + { + CoarseVector Ac(Coarse5d); + for (int k = 0; k < nbasis; k++) { + Aggregates.ProjectToSubspace(coarse_deflation_vecs[k], chi[k]); + RealD n = norm2(coarse_deflation_vecs[k]); + coarse_deflation_vecs[k] *= 1.0 / std::sqrt(n); + LinOpCoarse.HermOp(coarse_deflation_vecs[k], Ac); + coarse_deflation_evals[k] = real(TensorRemove(innerProduct(coarse_deflation_vecs[k], Ac))); + std::cout << GridLogMessage << "Coarse deflation eval[" << k << "] = " + << coarse_deflation_evals[k] << std::endl; + } + } + + // ── Coarse solve: CG + deflation guesser ────────────────────────────── + ConjugateGradient coarseCG(5.0e-2, 10000, false); + DeflatedGuesser coarseGuess(coarse_deflation_vecs, coarse_deflation_evals); + HPDSolver CoarseSolve(LinOpCoarse, coarseCG, coarseGuess); + + // ── ADEF2 outer solve ────────────────────────────────────────────────── + LatticeFermionD src(FGrid); random(RNG5, src); + LatticeFermionD result(FGrid); result = Zero(); + + TwoLevelADEF2 + HDCG(1.0e-8, 1000, + HermFineOp, + Smoother, + CoarseSolve, // used in PcgM1 + CoarseSolve, // used in Vstart + Aggregates); + + HDCG(src, result); + + // ── Reference RBCG ───────────────────────────────────────────────────── +#if 0 + { + SchurDiagMooeeOperator HermOpEO(Ddwf); + LatticeFermionD rb_src(FrbGrid); random(RNG5, rb_src); + LatticeFermionD rb_res(FrbGrid); rb_res = Zero(); + ConjugateGradient CG(1.0e-8, 30000, false); + CG(HermOpEO, rb_src, rb_res); + } +#endif + + Grid_finalize(); + return 0; +} diff --git a/examples/Example_pvdagm_mrhs.cc b/examples/Example_pvdagm_mrhs.cc new file mode 100644 index 000000000..cf22afa9c --- /dev/null +++ b/examples/Example_pvdagm_mrhs.cc @@ -0,0 +1,662 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./examples/Example_pvdagm_mrhs.cc + + 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. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution directory + *************************************************************************************/ + /* END LEGAL */ + +// MultiRHS (valence) two-level multigrid for PVdagM. +// +// Philosophy: NO block Krylov, NO per-RHS lockstep coefficients. A single +// GCR polynomial for the enlarged block-diagonal system diag(A,...,A): +// all inner products are summed over the RHS index, giving one alpha/beta +// per step shared by every RHS. +// +// Level structure mirrors Example_pvdagm: +// outer: MrhsPGCRNonHermitian on PVdagM over std::vector +// precon: V-cycle -- per-RHS fine post-smoother (16-step shifted GCR), +// batched restriction (MultiRHSBlockProject / GEMM), +// ONE coarse PGCR on the 6D mrhs coarse operator +// (MultiGeneralCoarsenedMatrix, GEMM mults -- the ~10x win), +// batched prolongation. +// +// The coarse operator is coarsened once with the standard single-RHS +// machinery (subspace cache reused) and imported via CopyMatrix. +// +// Memory note: outer restart history is 2*mmax*nrhs fine fields +// (49GB/field global at 48^3x96,Ls=24). Default mmax=8, nrhs=12 needs +// ~10TB for history; use OuterMmax / NRHS to fit the partition. +// +// Env vars: MASS, SUBSPACE_FILE, BLOCK (dotted e.g. 4.4.3.4), +// NRHS (default 12, multiple of Nsimd), +// FineSmootherShift, FineSmootherOrder, +// CoarseSolverTol, CoarseSolverOrder, +// OuterMmax, OuterNstep, OuterTol + +#include +#include +#include +#include + +using namespace std; +using namespace Grid; + +RealD FineSmootherShift = 0.1; +int FineSmootherOrder = 16; +RealD CoarseSolverTol = 0.03; +int CoarseSolverOrder = 200; +RealD OuterTol = 1.0e-8; +int OuterMmax = 8; +int OuterNstep = 8; +int Nrhs = 12; +RealD mass = 0.00078; + +void ParseEnvironment(void) +{ + if(getenv("MASS")) mass = atof(getenv("MASS")); + if(getenv("FineSmootherShift")) FineSmootherShift = atof(getenv("FineSmootherShift")); + if(getenv("FineSmootherOrder")) FineSmootherOrder = atoi(getenv("FineSmootherOrder")); + 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("NRHS")) Nrhs = atoi(getenv("NRHS")); + + std::cout << GridLogMessage << "PARAM: MASS " << mass << std::endl; + std::cout << GridLogMessage << "PARAM: NRHS " << Nrhs << std::endl; + std::cout << GridLogMessage << "PARAM: FineSmootherShift " << FineSmootherShift << std::endl; + std::cout << GridLogMessage << "PARAM: FineSmootherOrder " << FineSmootherOrder << std::endl; + std::cout << GridLogMessage << "PARAM: CoarseSolverTol " << CoarseSolverTol << std::endl; + std::cout << GridLogMessage << "PARAM: CoarseSolverOrder " << CoarseSolverOrder << std::endl; + std::cout << GridLogMessage << "PARAM: OuterTol " << OuterTol << std::endl; + std::cout << GridLogMessage << "PARAM: OuterMmax " << OuterMmax << std::endl; + std::cout << GridLogMessage << "PARAM: OuterNstep " << OuterNstep << std::endl; +} + +template +void saveSubspace(std::vector &subspace, std::string const fname){ +#ifdef HAVE_LIME + std::cout << Grid::GridLogMessage << "Saving subspace (" << subspace.size() << " vectors) to: " << fname << std::endl; + 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 + std::cout << Grid::GridLogMessage << "Loading subspace (" << subspace.size() << " vectors) from: " << fname << std::endl; + 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 +} + +template +class PVdagMLinearOperator : public LinearOperatorBase { + Matrix &_Mat; + Matrix &_PV; +public: + PVdagMLinearOperator(Matrix &Mat,Matrix &PV): _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); + } + 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 dot = innerProduct(in,out); + n1=real(dot); + n2=norm2(out); + } + void HermOp(const Field &in, Field &out){ + Field tmp(in.Grid()); + Op(in,tmp); + AdjOp(tmp,out); + } +}; + +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); + } +}; + +////////////////////////////////////////////////////////////////////// +// Minimal multi-RHS function interface (preconditioner slot) +////////////////////////////////////////////////////////////////////// +template +class MrhsLinearFunction { +public: + virtual void operator()(std::vector &in, std::vector &out) = 0; +}; + +////////////////////////////////////////////////////////////////////// +// Single-polynomial multi-RHS PGCR (non-Hermitian). +// +// Verbatim adaptation of PrecGeneralisedConjugateResidualNonHermitian +// to std::vector: every innerProduct / norm2 is SUMMED over the +// RHS index, so one alpha/beta per step is shared by all RHS -- the +// single GCR on the enlarged block-diagonal system. +////////////////////////////////////////////////////////////////////// +template +class MrhsPGCRNonHermitian { +public: + RealD Tolerance; + Integer MaxIterations; + int mmax; + int nstep; + int steps; + int level; + int ZeroGuess = 0; int FirstCycle = 0; // caller contract: zero guess => first-cycle r0 = src + 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; } + + /////////////////////////////////////////////////////////////// + // vector-of-fields linear algebra, reductions summed over rhs + /////////////////////////////////////////////////////////////// + static RealD vnorm2(std::vector &x){ + RealD s=0.0; for(auto &f : x) s+=norm2(f); return s; + } + static ComplexD vinnerProduct(std::vector &x, std::vector &y){ + ComplexD s(0.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 SolverTimer; + SolverTimer.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); + std::vector z (nrhs,grid); + std::vector Az(nrhs,grid); + + //////////////////////////////// + // history for flexible orthog: [mmax][nrhs] + //////////////////////////////// + std::vector< std::vector > q(mmax, std::vector(nrhs,grid)); + std::vector< std::vector > p(mmax, std::vector(nrhs,grid)); + std::vector qq(mmax); + + std::cout<(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); // never reached + return cp; + } +}; + +////////////////////////////////////////////////////////////////////// +// Trivial multi-RHS preconditioner +////////////////////////////////////////////////////////////////////// +template +class TrivialMrhsPrecon : public MrhsLinearFunction { +public: + virtual void operator()(std::vector &in, std::vector &out){ + for(int r=0;r<(int)in.size();r++) out[r]=in[r]; + } +}; + +////////////////////////////////////////////////////////////////////// +// MultiRHS two-level V-cycle. +// +// Mirrors MGPreconditioner in Example_pvdagm: +// out = in (trivial pre) [per rhs] +// r1 = in - A out [per rhs] +// batched blockProject -> pack -> ONE mrhs coarse PGCR -> unpack +// -> batched blockPromote; out += correction +// r2 = in - A out [per rhs] +// per-RHS fine post-smoother; out += smooth(r2) +////////////////////////////////////////////////////////////////////// +template +class MrhsTwoLevelMG : public MrhsLinearFunction { +public: + typedef MrhsCoarseVector CoarseVector; // same lattice type on Coarse5d and CoarseMrhs + + LinearOperatorBase &_FineOperator; + FineSmoother &_PostSmoother; // single-RHS smoother, looped + MultiRHSBlockProject &_Projector; + LinearFunction &_CoarseSolve; // PGCR on the 6D mrhs field + GridBase *_CoarseGrid; // Coarse5d (single rhs) + GridBase *_CoarseGridMrhs; // 6D + + 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(); + double t; + + std::vector vec1(nrhs,fgrid); + std::vector vec2(nrhs,fgrid); + + // Trivial pre-smoother: out = in (as in Example_pvdagm with simple_fine) + for(int r=0;rcoarse, pack rhs into 6D field + std::vector Csrc_split(nrhs,_CoarseGrid); + std::vector Csol_split(nrhs,_CoarseGrid); + CoarseVector CsrcMrhs(_CoarseGridMrhs); + CoarseVector CsolMrhs(_CoarseGridMrhs); + + t=-usecond(); + _Projector.blockProject(vec1,Csrc_split); + for(int r=0;rfine, add correction + t=-usecond(); + for(int r=0;r PVdagM_t; + typedef ShiftedPVdagMLinearOperator ShiftedPVdagM_t; + typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef MultiGeneralCoarsenedMatrix MrhsLittleDiracOperator; + typedef LittleDiracOperator::CoarseVector CoarseVector; + typedef Aggregation Subspace; + + PVdagM_t PVdagM(Ddwf,Dpv); + ShiftedPVdagM_t ShiftedPVdagM(FineSmootherShift,Ddwf,Dpv); + + NextToNearestStencilGeometry5D geom(Coarse5d); + + //////////////////////////////////////////////////////////// + // Subspace: load from cache or generate + //////////////////////////////////////////////////////////// + std::string subspace_file = "/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle/subspace_nb" + + std::to_string(nbasis) + ".scidac"; + if ( getenv("SUBSPACE_FILE") ) subspace_file = std::string(getenv("SUBSPACE_FILE")); + + uint64_t file_exists = 0; + if ( UGrid->IsBoss() ) { std::ifstream f(subspace_file); file_exists = f.good() ? 1 : 0; } + UGrid->GlobalSum(file_exists); + + const int cb = 0; + Subspace AggregatesGCR(Coarse5d,FGrid,cb); + + if ( file_exists ) { + std::cout << GridLogMessage << "*** Loading subspace from disk ***" << std::endl; + loadSubspace(AggregatesGCR.subspace, subspace_file); + } else { + std::cout << GridLogMessage << "*** GCR subspace generation ***" << std::endl; + AggregatesGCR.CreateSubspaceGCR(RNG5,PVdagM,nbasis); + saveSubspace(AggregatesGCR.subspace, subspace_file); + } + + //////////////////////////////////////////////////////////// + // Coarsen once (single-RHS machinery), import into mrhs op. + // NB CoarsenOperator block-orthogonalises subspace in place; + // ImportBasis AFTER so the projector matches the coarse op. + //////////////////////////////////////////////////////////// + MrhsLittleDiracOperator mrhsLittleDiracOpPV(geom,CoarseMrhs); + { + // Scope the single-RHS operator so its padded _A is freed after import. + // At small local volumes the depth-2 padded cell inflates ~8x + // (e.g. 2^4 blocking on 432 ranks: local 8x4x4x12 -> padded 12x8x8x16, + // ~23GB/GCD for _A alone -> OOM if kept alive). + LittleDiracOperator LittleDiracOpPV(geom,FGrid,Coarse5d); + LittleDiracOpPV.CoarsenOperator(PVdagM, AggregatesGCR); + mrhsLittleDiracOpPV.CopyMatrix(LittleDiracOpPV); + } + + MultiRHSBlockProject MrhsProjector; + MrhsProjector.Allocate(nbasis,FGrid,Coarse5d); + MrhsProjector.ImportBasis(AggregatesGCR.subspace); + + //////////////////////////////////////////////////////////// + // Solvers + //////////////////////////////////////////////////////////// + NonHermitianLinearOperator mrhsLinOpCoarse(mrhsLittleDiracOpPV); + + TrivialPrecon simpleC; + PrecGeneralisedConjugateResidualNonHermitian + L2PGCRmrhs(CoarseSolverTol, CoarseSolverOrder/20, mrhsLinOpCoarse, simpleC, 20, 20); + L2PGCRmrhs.Level(2); + L2PGCRmrhs.Name("Couter"); + L2PGCRmrhs.SetZeroGuess(1); // caller zeroes CsolMrhs + + TrivialPrecon simple_fine; + PrecGeneralisedConjugateResidualNonHermitian + SmootherGCR(0.0, 1, ShiftedPVdagM, simple_fine, FineSmootherOrder, FineSmootherOrder); + SmootherGCR.Level(1); + SmootherGCR.Name("Fsmoother"); + SmootherGCR.SetZeroGuess(1); // caller zeroes vec2[r] + + typedef PrecGeneralisedConjugateResidualNonHermitian FineSmoother_t; + MrhsTwoLevelMG + TwoLevelPrecon(PVdagM, SmootherGCR, MrhsProjector, L2PGCRmrhs, Coarse5d, CoarseMrhs); + + MrhsPGCRNonHermitian + L1PGCRmrhs(OuterTol, 1000, PVdagM, TwoLevelPrecon, OuterMmax, OuterNstep); + L1PGCRmrhs.Level(1); + L1PGCRmrhs.Name("Fouter"); + L1PGCRmrhs.SetZeroGuess(1); // sol[r]=Zero() at source setup + + //////////////////////////////////////////////////////////// + // Sources and solve + //////////////////////////////////////////////////////////// + std::vector src(nrhs,FGrid); + std::vector sol(nrhs,FGrid); + 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. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ + +// +// 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. +// +// Differences from Example_pvdagm_mrhs_3level_DenseCoarseMatrix.cc: +// +// * The coarse space is UNVECTORISED (sComplexD). The fine space stays +// vectorised. MultiRHSBlockProject carries the mixed layout. +// +// * One operator, not two. V1 needed GeneralCoarsenedMatrix to coarsen and +// MultiGeneralCoarsenedMatrix to apply, bridged by CopyMatrix. V2 does +// both, and single versus multiRHS is SetGrid on the same object with the +// matrix elements built once. +// +// * Nrhs is unconstrained. V1 required nrhs % vComplex::Nsimd() == 0 because +// its multiRHS grid carried the SIMD in the rhs direction. +// +// * CoarsenOperator takes the subspace vectors, not an Aggregation. It block +// orthonormalises them IN PLACE -- the vectors are far too large to copy +// defensively -- so rawNull is taken first and the RAW vectors are what +// define the L2 null space. Do not insert an Orthogonalise() anywhere: +// projecting a block-orthonormal vector onto its own block-orthonormalised +// 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 +// + +#include +#include +#include +#include + +#include + +using namespace std; +using namespace Grid; + +// Compile time so it can be cut down for laptop runs: -DNBASIS=8 +#ifndef NBASIS +#define NBASIS 60 +#endif + +RealD mass = 0.00078; +int Nrhs = 12; +int Ls = 24; +int CoarsenBatch = 9; +std::vector lat_size({48,48,48,96}); + +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("LATT")){ + Coordinate l; + GridCmdOptionIntVector(std::string(getenv("LATT")),l); + GRID_ASSERT(l.size()==4); + for(int d=0;d<4;d++) lat_size[d]=l[d]; + } + + std::cout << GridLogMessage << "PARAM: LATT " + << lat_size[0]<<"."< +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 +} + +////////////////////////////////////////////////////////////////////// +// A = PV^dag M (non-Hermitian) +////////////////////////////////////////////////////////////////////// +template +class PVdagMLinearOperator : public LinearOperatorBase { + Matrix &_Mat; Matrix &_PV; +public: + PVdagMLinearOperator(Matrix &Mat,Matrix &PV): _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); } + 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); } +}; + +////////////////////////////////////////////////////////////////////// +// || - 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. +////////////////////////////////////////////////////////////////////// +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); +} + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + ParseEnvironment(); + + RealD M5=1.8, b=1.5, c=0.5; + const int nbasis=NBASIS; + const int nrhs=Nrhs; + const int batch=CoarsenBatch; + + Coordinate mpi = GridDefaultMpi(); + Coordinate fsimd= GridDefaultSimd(Nd,vComplex::Nsimd()); + + GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(lat_size,fsimd,mpi); + GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); + GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); + + // Level 1 blocking (default 2^4) + Coordinate clatt = lat_size; + Coordinate Block({2,2,2,2}); + if ( getenv("BLOCK") ){ GridCmdOptionIntVector(std::string(getenv("BLOCK")),Block); GRID_ASSERT(Block.size()==4); } + for(int d=0;d<4;d++){ GRID_ASSERT(lat_size[d]%Block[d]==0); clatt[d]=lat_size[d]/Block[d]; } + std::cout << GridLogMessage << "Block " << Block << " coarse lattice " << clatt << std::endl; + + ////////////////////////////////////////////////////////////////////// + // The coarse space is unvectorised. The 5D coarse grid is built here + // rather than through SpaceTimeGrid so the SIMD layout is ours. + ////////////////////////////////////////////////////////////////////// + Coordinate c5latt({1,clatt[0],clatt[1],clatt[2],clatt[3]}); + Coordinate c5simd({1,1,1,1,1}); + Coordinate c5mpi ({1,mpi[0],mpi[1],mpi[2],mpi[3]}); + GridCartesian *Coarse5d = new GridCartesian(c5latt,c5simd,c5mpi); + + // 6D coarse multiRHS grid: rhs is dim 0, undistributed and unvectorised. + // No divisibility constraint on nrhs, unlike V1. + Coordinate cmlatt({nrhs,1,clatt[0],clatt[1],clatt[2],clatt[3]}); + Coordinate cmsimd({1,1,1,1,1,1}); + Coordinate cmmpi ({1,1,mpi[0],mpi[1],mpi[2],mpi[3]}); + GridCartesian *CoarseMrhs = new GridCartesian(cmlatt,cmsimd,cmmpi); + + // 6D coarse grid at the coarsening batch, used only while CoarsenOperator + // runs. The matrix elements survive the change back to nrhs. + Coordinate cblatt({batch,1,clatt[0],clatt[1],clatt[2],clatt[3]}); + GridCartesian *CoarseBatch = new GridCartesian(cblatt,cmsimd,cmmpi); + + // 6D fine grid carrying the coarsening batch: fine SIMD layout preserved + Coordinate fmlatt({batch,Ls,lat_size[0],lat_size[1],lat_size[2],lat_size[3]}); + Coordinate fmsimd({1,1,fsimd[0],fsimd[1],fsimd[2],fsimd[3]}); + Coordinate fmmpi ({1,1,mpi[0],mpi[1],mpi[2],mpi[3]}); + GridCartesian *FineMrhs = new GridCartesian(fmlatt,fmsimd,fmmpi); + + std::cout << GridLogMessage << "Nsimd fine " << FGrid->Nsimd() + << " coarse " << Coarse5d->Nsimd() << std::endl; + + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers({1,2,3,4}); + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers({5,6,7,8}); + + ////////////////////////////////////////////////////////////////////// + // Gauge field + ////////////////////////////////////////////////////////////////////// + LatticeGaugeField Umu(UGrid); + if ( getenv("HOT_START") ) { + std::cout << GridLogMessage << "Hot start gauge field" << std::endl; + SU::HotConfiguration(RNG4,Umu); + } else { + std::string file("/ccs/home/poare/ckpoint_lat.1000"); + if ( getenv("CONFIG") ) file = std::string(getenv("CONFIG")); + std::cout << GridLogMessage << "Reading gauge field " << file << std::endl; + FieldMetaData header; + NerscIO::readConfiguration(Umu,header,file); + } + + MobiusFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c); + MobiusFermionD Dpv (Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,1.0, M5,b,c); + + typedef PVdagMLinearOperator PVdagM_t; + PVdagM_t PVdagM(Ddwf,Dpv); + + ////////////////////////////////////////////////////////////////////// + // Level 1 types: unvectorised coarse scalar + ////////////////////////////////////////////////////////////////////// + typedef sTComplexD CComplexS; + typedef MultiGeneralCoarsenedOperatorV2 CoarseOperator; + typedef CoarseOperator::CoarseVector CoarseVector; + typedef Aggregation Subspace; + + NextToNearestStencilGeometry5D geom(Coarse5d); + + ////////////////////////////////////////////////////////////////////// + // Subspace: load RAW (no Orthogonalise!), or generate. + // + // The Aggregation is scaffolding for CreateSubspaceGCR only. That runs + // entirely on the fine grid and ends in GlobalOrthonormalise, which is a + // whole-lattice Gram-Schmidt, so the coarse grid it holds is never + // dereferenced and may be the unvectorised one. + ////////////////////////////////////////////////////////////////////// + std::string subspace_file = "subspace_nb" + std::to_string(nbasis) + ".scidac"; + if ( getenv("SUBSPACE_FILE") ) subspace_file = std::string(getenv("SUBSPACE_FILE")); + uint64_t file_exists=0; + if ( UGrid->IsBoss() ){ std::ifstream f(subspace_file); file_exists=f.good()?1:0; } + UGrid->GlobalSum(file_exists); + + const int cb=0; + Subspace AggregatesGCR(Coarse5d,FGrid,cb); + if ( file_exists ){ + std::cout << GridLogMessage << "*** Loading subspace from disk (kept RAW) ***" << std::endl; + loadSubspace(AggregatesGCR.subspace, subspace_file); + } else { + std::cout << GridLogMessage << "*** GCR subspace generation ***" << std::endl; + AggregatesGCR.CreateSubspaceGCR(RNG5,PVdagM,nbasis); + saveSubspace(AggregatesGCR.subspace, subspace_file); + } + + // RAW copy BEFORE CoarsenOperator block-orthonormalises in place. + std::vector rawNull(nbasis,FGrid); + for(int k=0;k MrhsPVdagM(PVdagM,FGrid,batch); + CoarseOpPV.CoarsenOperator(MrhsPVdagM,FineMrhs,AggregatesGCR.subspace,Coarse5d); + } else { + // PVdagM is single RHS: apply it directly, batch on the coarse side. + CoarseOpPV.CoarsenOperator(PVdagM,AggregatesGCR.subspace,Coarse5d,batch); + } + + // Stay on the batch grid: the L2 coarsening drives this operator at the + // batch. It is switched to the solve Nrhs once L2 is built. + + ////////////////////////////////////////////////////////////////////// + // psi_coarse = P^dag (RAW fine null) -> Galerkin images that carry the + // near null content, and are free: A_c (P psi) = P A psi. + ////////////////////////////////////////////////////////////////////// + MultiRHSBlockProject MrhsProjector; + MrhsProjector.Allocate(nbasis,FGrid,Coarse5d); + MrhsProjector.ImportBasis(AggregatesGCR.subspace); // block orthonormal basis + + std::vector psi_coarse(nbasis,Coarse5d); + 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 ); + } + + ////////////////////////////////////////////////////////////////////// + // Optional cross check of the coarse matrix elements against the V1 + // path, which needs a vectorised coarse space. Block Gram-Schmidt is + // idempotent, so V1 may re-orthonormalise the same vectors in place + // without a second copy of the subspace. + ////////////////////////////////////////////////////////////////////// + if ( getenv("V1_CHECK") ) { + + typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef MultiGeneralCoarsenedMatrix MrhsLittleDiracOperator; + typedef Aggregation SubspaceV; + + Coordinate v5latt({1,clatt[0],clatt[1],clatt[2],clatt[3]}); + Coordinate v5simd({1,fsimd[0],fsimd[1],fsimd[2],fsimd[3]}); + GridCartesian *Coarse5dV = new GridCartesian(v5latt,v5simd,c5mpi); + + int nrhs_v1 = vComplex::Nsimd(); + Coordinate vmlatt({nrhs_v1,1,clatt[0],clatt[1],clatt[2],clatt[3]}); + Coordinate vmsimd({vComplex::Nsimd(),1,1,1,1,1}); + GridCartesian *CoarseMrhsV = new GridCartesian(vmlatt,vmsimd,cmmpi); + + NextToNearestStencilGeometry5D geomV(Coarse5dV); + + SubspaceV AggV(Coarse5dV,FGrid,cb); + for(int k=0;k h1(sites),h2(sites); + acceleratorCopyFromDevice(&mrhsV1.BLAS_A[p][0], &h1[0],sites*sizeof(calcMatrix)); + acceleratorCopyFromDevice(&CoarseOpPV.BLAS_A[p][0],&h2[0],sites*sizeof(calcMatrix)); + ComplexD *w1=(ComplexD *)&h1[0]; + ComplexD *w2=(ComplexD *)&h2[0]; + int64_t words = sites*sizeof(calcMatrix)/sizeof(ComplexD); + for(int64_t i=0;i CComplexS2; + typedef MultiGeneralCoarsenedOperatorV2 CoarseCoarseOperator; + typedef CoarseCoarseOperator::CoarseVector CoarseCoarseVector; + + NextToNearestStencilGeometry5D geom2(CoarseCoarse5d); + + // RAW copy of the coarse null vectors, for the same reason as rawNull: + // the L2 CoarsenOperator block-orthonormalises its subspace in place, and + // the L3 basis must be defined by the vectors that still carry content. + std::vector rawPsi(nbasis,Coarse5d); + for(int k=0;k LinOpCoarse(CoarseOpPV); + + std::cout << GridLogMessage << "*** L2 CoarsenOperator, batch "< MrhsProjectorL2; + MrhsProjectorL2.Allocate(nbasis,Coarse5d,CoarseCoarse5d); + MrhsProjectorL2.ImportBasis(psi_coarse); // block orthonormal basis + + { + std::vector psi_cc(nbasis,CoarseCoarse5d); + MrhsProjectorL2.blockProject(rawPsi,psi_cc); // RAW vectors in + + RealD defect = GramDefect(psi_cc); + RealD leak = std::sqrt((double)CoarseCoarse5d->gSites()); + std::cout << GridLogMessage << "GUARD: || - I||_F = " << defect + << " (~0.23 good; ~sqrt(N_cc)=" << leak << " = e_k leak)" << std::endl; + GRID_ASSERT( defect < leak ); + } + rawPsi.clear(); rawPsi.shrink_to_fit(); + + ////////////////////////////////////////////////////////////////////// + // Both coarse operators apply on their solve grids + ////////////////////////////////////////////////////////////////////// + { + GridParallelRNG cRNG(Coarse5d); cRNG.SeedFixedIntegers({3,4,5,6}); + CoarseVector cin(CoarseMrhs), cout_(CoarseMrhs); + random(cRNG,cin); + CoarseOpPV.M(cin,cout_); + std::cout << GridLogMessage << "L1 apply |in|^2 = " << norm2(cin) + << " |M in|^2 = " << norm2(cout_) << std::endl; + GRID_ASSERT( norm2(cout_) > 0.0 ); + + GridParallelRNG ccRNG(CoarseCoarse5d); ccRNG.SeedFixedIntegers({7,8,9,10}); + CoarseCoarseVector ccin(CoarseCoarseMrhs), ccout(CoarseCoarseMrhs); + random(ccRNG,ccin); + CoarseOpL2.M(ccin,ccout); + std::cout << GridLogMessage << "L2 apply |in|^2 = " << norm2(ccin) + << " |M in|^2 = " << norm2(ccout) << std::endl; + GRID_ASSERT( norm2(ccout) > 0.0 ); + } + + std::cout << GridLogMessage << "*** stage two complete: L1 and L2 coarse operators built ***" << std::endl; + + Grid_finalize(); +}