diff --git a/examples/Example_pvdagm_3level.cc b/examples/Example_pvdagm_3level.cc new file mode 100644 index 000000000..f7c4387b6 --- /dev/null +++ b/examples/Example_pvdagm_3level.cc @@ -0,0 +1,699 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./examples/Example_pvdagm_3level.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 + +#include +#include +#include + +using namespace std; +using namespace Grid; + +template void readFile(T& out, std::string const fname){ + #ifdef HAVE_LIME + std::cout << Grid::GridLogMessage << "Reading: " << fname << std::endl; + Grid::emptyUserRecord record; + Grid::ScidacReader SR; + SR.open(fname); + SR.readScidacFieldRecord(out, record); + SR.close(); + #endif +} + +template void writeFile(T& in, std::string const fname){ + #ifdef HAVE_LIME + std::cout << Grid::GridLogMessage << "Writing: " << fname << std::endl; + Grid::emptyUserRecord record; + Grid::ScidacWriter SW(in.Grid()->IsBoss()); + SW.open(fname); + SW.writeScidacFieldRecord(in, record); + SW.close(); + #endif +} + +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; + int nApp; + int nAppDag; +public: + PVdagMLinearOperator(Matrix &Mat,Matrix &PV): _Mat(Mat),_PV(PV), nApp(0), nAppDag(0) {}; + + 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); + nApp++; + } + void AdjOp (const Field &in, Field &out){ + Field tmp(in.Grid()); + _PV.M(in,tmp); + _Mat.Mdag(tmp,out); + nAppDag++; + } + void clear() { nApp = 0; nAppDag = 0; } + void getApplications() { + std::cout << GridLogMessage << "# applications of PVdagM: " << nApp << std::endl; + std::cout << GridLogMessage << "# applications of PVdagM^dag: " << nAppDag << std::endl; + std::cout << GridLogMessage << "# applications total: " << nApp + nAppDag << std::endl; + } + 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 MdagPVLinearOperator : public LinearOperatorBase { + Matrix &_Mat; + Matrix &_PV; +public: + MdagPVLinearOperator(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()); + _PV.M(in,tmp); + _Mat.Mdag(tmp,out); + } + void AdjOp (const Field &in, Field &out){ + Field tmp(in.Grid()); + _Mat.M(in,tmp); + _PV.Mdag(tmp,out); + } + void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ + 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; + RealD shift; +public: + 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); + } +}; + +// Lüscher deflated guesser (arXiv:0706.2298 Sec A.3) for a non-Hermitian solve. +// C_{st} = ; guess = sum_s c_s psi[s] where c = C^{-1} psi† src. +template +class LuscherGuesser : public LinearFunction { + const std::vector ψ + Eigen::MatrixXcd C_inv; +public: + using LinearFunction::operator(); + LuscherGuesser(const std::vector &psi_, const Eigen::MatrixXcd &Cinv_) + : psi(psi_), C_inv(Cinv_) {} + virtual void operator()(const Field &src, Field &guess) { + int N = psi.size(); + Eigen::VectorXcd b(N); + for (int t = 0; t < N; t++) + b(t) = TensorRemove(innerProduct(psi[t], src)); + Eigen::VectorXcd c = C_inv * b; + guess = Zero(); + for (int s = 0; s < N; s++) + guess += ComplexD(c(s)) * psi[s]; + } +}; + +template +class MGPreconditioner : public LinearFunction< Lattice > { +public: + using LinearFunction >::operator(); + + typedef Aggregation Aggregates; + typedef typename Aggregation::FineField FineField; + typedef typename Aggregation::CoarseVector CoarseVector; + typedef typename Aggregation::CoarseMatrix CoarseMatrix; + typedef LinearOperatorBase FineOperator; + typedef LinearFunction FineSmoother; + typedef LinearOperatorBase CoarseOperator; + typedef LinearFunction CoarseSolver; + + Aggregates & _Aggregates; + FineOperator & _FineOperator; + FineSmoother & _PreSmoother; + FineSmoother & _PostSmoother; + CoarseOperator & _CoarseOperator; + CoarseSolver & _CoarseSolve; + CoarseSolver & _CoarseGuesser; + + int level; void Level(int lv) {level = lv; }; + + MGPreconditioner(Aggregates &Agg, + FineOperator &Fine, + FineSmoother &PreSmoother, + FineSmoother &PostSmoother, + CoarseOperator &CoarseOperator_, + CoarseSolver &CoarseSolve_, + CoarseSolver &CoarseGuesser_) + : _Aggregates(Agg), + _FineOperator(Fine), + _PreSmoother(PreSmoother), + _PostSmoother(PostSmoother), + _CoarseOperator(CoarseOperator_), + _CoarseSolve(CoarseSolve_), + _CoarseGuesser(CoarseGuesser_), + level(1) { } + + virtual void operator()(const FineField &in, FineField & out) + { + GridBase *CoarseGrid = _Aggregates.CoarseGrid; + CoarseVector Csrc(CoarseGrid); + CoarseVector Csol(CoarseGrid); + FineField vec1(in.Grid()); + FineField vec2(in.Grid()); + + double t; + out = Zero(); + t=-usecond(); + _PreSmoother(in,out); + t+=usecond(); + std::cout< +class ShiftedLinearOperator : public LinearOperatorBase { + LinearOperatorBase &_Op; + RealD shift; +public: + ShiftedLinearOperator(RealD _shift, LinearOperatorBase &Op) : shift(_shift), _Op(Op) {} + 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); } +}; + +template +void runMG( + GridCartesian *FGrid, + GridCartesian *Coarse5d, + GridCartesian *CoarseCoarse5d, + NextToNearestStencilGeometry5D geom, + PVdagM_t &PVdagM, + ShiftedPVdagM_t &ShiftedPVdagM, + Subspace &AggregatesPD +) { + std::vector subspace = AggregatesPD.subspace; + assert((int)subspace.size() == NB); + const int nbasis = NB; + const int cb = 0; + + CoarseVector c_src(Coarse5d); + CoarseVector c_res(Coarse5d); + Complex one(1.0); + + LatticeFermionD f_src(FGrid); + LatticeFermionD f_res(FGrid); + + TrivialPrecon simpleC; + TrivialPrecon simple_fine; + + ////////////////////////////////////////////////////////////////////// + // Level 0→1: coarsen PVdagM, build LinOpCoarse + ////////////////////////////////////////////////////////////////////// + LittleDiracOperator LittleDiracOpPV(geom, FGrid, Coarse5d); + LittleDiracOpPV.CoarsenOperator(PVdagM, AggregatesPD); + + NonHermitianLinearOperator LinOpCoarse(LittleDiracOpPV); + + ////////////////////////////////////////////////////////////////////// + // Baseline: plain PGCR on LinOpCoarse (reference for comparison) + ////////////////////////////////////////////////////////////////////// + std::cout< L2PGCR_baseline(3.0e-2,1100,LinOpCoarse,simpleC,10,10); + L2PGCR_baseline.Level(2); + L2PGCR_baseline.Name("Cbaseline"); + c_src = one; + c_res = Zero(); + L2PGCR_baseline(c_src,c_res); + + ////////////////////////////////////////////////////////////////////// + // psi_coarse: coarse projections of pre-GS fine null vectors. + // These are the Level 1 near-null vectors, promoted from Level 0. + // Used as the aggregation basis for Level 1→2 coarsening. + ////////////////////////////////////////////////////////////////////// + std::vector psi_coarse(nbasis, Coarse5d); + for (int k = 0; k < nbasis; k++) + AggregatesPD.ProjectToSubspace(psi_coarse[k], subspace[k]); + + ////////////////////////////////////////////////////////////////////// + // Diagnostics: W (fine projected matrix) and C (Galerkin check) + ////////////////////////////////////////////////////////////////////// + { + Eigen::MatrixXcd W = Eigen::MatrixXcd::Zero(nbasis, nbasis); + LatticeFermion ftmp(FGrid); + for (int j = 0; j < nbasis; j++) { + PVdagM.Op(subspace[j], ftmp); + for (int i = 0; i < nbasis; i++) + W(i,j) = TensorRemove(innerProduct(subspace[i], ftmp)); + } + RealD normW = W.norm(); + std::cout << GridLogMessage << "Fine projected matrix ||W|| = " << normW << std::endl; + + Eigen::MatrixXcd C = Eigen::MatrixXcd::Zero(nbasis, nbasis); + CoarseVector Ac(Coarse5d); + for (int l = 0; l < nbasis; l++) { + LinOpCoarse.Op(psi_coarse[l], Ac); + for (int k = 0; k < nbasis; k++) + C(k,l) = TensorRemove(innerProduct(psi_coarse[k], Ac)); + } + RealD normC = C.norm(); + RealD normCmCdag = (C - C.adjoint()).norm(); + std::cout << GridLogMessage << "Coarse null matrix ||C|| = " << normC << std::endl; + std::cout << GridLogMessage << "Coarse null matrix ||C - C†||/||C|| = " << normCmCdag/normC << std::endl; + std::cout << GridLogMessage << "Galerkin check ||C||/||W|| = " << normC/normW << std::endl; + } + + ////////////////////////////////////////////////////////////////////// + // Level 1→2: set up aggregation using psi_coarse as subspace. + // Block factor 2,2,3,2 (removes odd local sublattice in z given MPI + // geometry 3×6×4×4 where z-local at Level 1 is 6). + // psi_coarse are assigned directly; CoarsenOperator performs + // block-GS orthogonalisation before building LinOpCoarseCoarse. + ////////////////////////////////////////////////////////////////////// + // innerProduct(CoarseSiteObj, CoarseSiteObj) returns iScalar, so CComplex + // for the L1→L2 level must be iScalar, not vTComplex. + typedef typename CoarseVector::vector_object CoarseSiteObj; + typedef iScalar vTTComplex; + typedef GeneralCoarsenedMatrix LittleDiracOperatorL2; + typedef typename LittleDiracOperatorL2::CoarseVector CoarseCoarseVector; + typedef Aggregation SubspaceL2; + typedef MGPreconditioner L1to2MG; + + SubspaceL2 AggregatesL2(CoarseCoarse5d, Coarse5d, cb); + for (int k = 0; k < nbasis; k++) + AggregatesL2.subspace[k] = psi_coarse[k]; + + NextToNearestStencilGeometry5D geom2(CoarseCoarse5d); + LittleDiracOperatorL2 LittleDiracOpL2(geom2, Coarse5d, CoarseCoarse5d); + LittleDiracOpL2.CoarsenOperator(LinOpCoarse, AggregatesL2); + + NonHermitianLinearOperator LinOpCC(LittleDiracOpL2); + + TrivialPrecon simpleCC; + + ////////////////////////////////////////////////////////////////////// + // Lüscher deflation guesser for L3PGCR. + // Step 1: project psi_coarse[k] (promoted fine null vectors) to + // CoarseCoarseVector space — these cover the zero-momentum + // component of the near-null space of LinOpCC. + // Step 2: breed Nextra additional null vectors directly on LinOpCC + // using GCR with random sources — these pick up near-null + // modes at all spatial frequencies not spanned by step 1. + // Step 3: build C_{st} = over the + // full augmented basis and invert directly via Eigen LU. + ////////////////////////////////////////////////////////////////////// + std::vector psi_cc(nbasis, CoarseCoarse5d); + for (int k = 0; k < nbasis; k++) + AggregatesL2.ProjectToSubspace(psi_cc[k], psi_coarse[k]); + + { + int Nextra = nbasis; // breed as many extra as we have promoted ones + if ( getenv("CC_NEXTRA") ) Nextra = atoi(getenv("CC_NEXTRA")); + GridParallelRNG RNG_CC(CoarseCoarse5d); + RNG_CC.SeedFixedIntegers({11,13,17,19}); + PrecGeneralisedConjugateResidualNonHermitian + nullGCR(1e-2, 200, LinOpCC, simpleCC, 32, 32); + CoarseCoarseVector tmp(CoarseCoarse5d); + for (int k = 0; k < Nextra; k++) { + CoarseCoarseVector src(CoarseCoarse5d); + gaussian(RNG_CC, src); + tmp = Zero(); + nullGCR(src, tmp); + psi_cc.push_back(tmp); + } + std::cout << GridLogMessage << "LinOpCC deflation basis: " << nbasis + << " promoted + " << Nextra << " bred = " << psi_cc.size() << " total" << std::endl; + } + + const int Naug = psi_cc.size(); + Eigen::MatrixXcd Ccc = Eigen::MatrixXcd::Zero(Naug, Naug); + { + CoarseCoarseVector Acc(CoarseCoarse5d); + for (int l = 0; l < Naug; l++) { + LinOpCC.Op(psi_cc[l], Acc); + for (int k = 0; k < Naug; k++) + Ccc(k,l) = TensorRemove(innerProduct(psi_cc[k], Acc)); + } + } + { + RealD normCcc = Ccc.norm(); + RealD normCccmCdag = (Ccc - Ccc.adjoint()).norm(); + std::cout << GridLogMessage << "Coarse-coarse deflation matrix ||Ccc|| = " << normCcc << std::endl; + std::cout << GridLogMessage << "Coarse-coarse deflation matrix ||Ccc-Ccc†||/||Ccc|| = " << normCccmCdag/normCcc << std::endl; + } + Eigen::MatrixXcd Ccc_inv = Ccc.inverse(); + LuscherGuesser CCDeflGuesser(psi_cc, Ccc_inv); + + ////////////////////////////////////////////////////////////////////// + // Level 2 solver: plain GCR, no further coarsening + ////////////////////////////////////////////////////////////////////// + PrecGeneralisedConjugateResidualNonHermitian L3PGCR(1.0e-1,200,LinOpCC,simpleCC,16,16); + L3PGCR.Level(3); + L3PGCR.Name("CCouter"); + + ////////////////////////////////////////////////////////////////////// + // Coarse-level GCR smoother for Level 1→2 V-cycle. + // Mirrors fine-grid SmootherGCR: shifted operator + fixed step count. + // coarse_smoother_shift and coarse_smoother_nstep are the tuning knobs. + ////////////////////////////////////////////////////////////////////// + RealD coarse_smoother_shift = 0.01; + int coarse_smoother_nstep = 12; + if(getenv("coarse_smoother_shift")) coarse_smoother_shift = atof(getenv("coarse_smoother_shift")); + if(getenv("coarse_smoother_nstep")) coarse_smoother_nstep = atoi(getenv("coarse_smoother_nstep")); + + ShiftedLinearOperator ShiftedLinOpCoarse(coarse_smoother_shift, LinOpCoarse); + PrecGeneralisedConjugateResidualNonHermitian CoarseSmootherGCR(0.01,1,ShiftedLinOpCoarse,simpleC,coarse_smoother_nstep,coarse_smoother_nstep); + CoarseSmootherGCR.SetZeroGuess(1); // smoother slot: caller zeroes guess + CoarseSmootherGCR.Level(2); + CoarseSmootherGCR.Name("Csmoother"); + + ////////////////////////////////////////////////////////////////////// + // Level 1→2 V-cycle preconditioner. + ////////////////////////////////////////////////////////////////////// + L1to2MG L1to2Precon(AggregatesL2, + LinOpCoarse, + simpleC, // no pre-smoother (matches fine-grid setup) + CoarseSmootherGCR, // post-smoother: 12 GCR steps + LinOpCC, + L3PGCR, + CCDeflGuesser); // Lüscher guesser: psi_cc C^{-1} psi_cc† + + ////////////////////////////////////////////////////////////////////// + // Standalone Level 1 two-level solve test. + // Compare against plain PGCR baseline above. + ////////////////////////////////////////////////////////////////////// + std::cout< L2MGsolver(3.0e-2,200,LinOpCoarse,L1to2Precon,16,16); + L2MGsolver.Level(2); + L2MGsolver.Name("Couter"); + c_res = Zero(); + L2MGsolver(c_src,c_res); + + std::cout << GridLogMessage << "Level 1 two-level test: PVdagM operator uses:" << std::endl; + PVdagM.getApplications(); + PVdagM.clear(); + + ////////////////////////////////////////////////////////////////////// + // Full three-level outer solve + ////////////////////////////////////////////////////////////////////// + std::cout< SmootherGCR(0.01,1,ShiftedPVdagM,simple_fine,16,16); + SmootherGCR.SetZeroGuess(1); // pre/post smoother slots zero their guess + SmootherGCR.Level(1); + SmootherGCR.Name("Fsmoother"); + + f_src = one; + + // Pre-smoother: none (TrivialPrecon); post-smoother: shifted PGCR. + // Coarse solver: L2MGsolver (PGCR preconditioned by Level 1→2 V-cycle). + TwoLevelMG ThreeLevelPrecon(AggregatesPD, + PVdagM, + simple_fine, + SmootherGCR, + LinOpCoarse, + L2MGsolver, + simpleC); + + PrecGeneralisedConjugateResidualNonHermitian L1PGCR(1.0e-8,1000,PVdagM,ThreeLevelPrecon,16,16); + L1PGCR.Level(1); + L1PGCR.Name("Fouter"); + + f_res = Zero(); + L1PGCR(f_src,f_res); + + std::cout << GridLogMessage << "Three-level outer solve: PVdagM operator uses:" << std::endl; + PVdagM.getApplications(); + PVdagM.clear(); +} + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + const int Ls = 24; + RealD M5 = 1.8; + RealD b = 1.5; + RealD c = 0.5; + RealD mass = 0.00078; + if ( getenv("MASS") ) mass = atof(getenv("MASS")); + + const int nbasis = 60; + + std::cout << GridLogMessage << "Mass: " << mass << ", Ls: " << Ls << ", b=" << b << ", c=" << c << std::endl; + std::cout << GridLogMessage << "nbasis: " << nbasis << std::endl; + + 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); + + // Level 1 coarse grid: block 2^4 from fine (48×48×48×96 → 24×24×24×48, Ls=1) + Coordinate clatt = lat_size; + for (int d = 0; d < 4; d++) clatt[d] /= 2; + std::cout << GridLogMessage << "Level 1 coarse lattice: " << clatt << std::endl; + + GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); + GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d); + + // Level 2 coarse-coarse grid: block 2,2,3,3 from Level 1 (24×24×24×48 → 12×12×8×16, Ls=1). + // MPI geometry 3.6.4.4 (288 ranks): fine local {16,8,12,24}. + // Level 1 local {8,4,6,12}; Level 2 local {4,2,2,4}. + // z blocked by 3: z-Level1-local=6; 6/3=2 (even), 6/2=3 (odd) → must use 3. + // t blocked by 3: t-Level1-local=12; 12/3=4 divisible by Nsimd=4 (gen-simd-width=64). + // t-block=2 gives t2-local=6, 6 mod 4 ≠ 0, fails Grid SIMD assertion. ✓ + // With {4,2,2,4}: Nsimd=4 goes into x or t (both =4). ✓ + Coordinate clatt2 = clatt; + clatt2[0] /= 2; + clatt2[1] /= 2; + clatt2[2] /= 3; + clatt2[3] /= 3; + std::cout << GridLogMessage << "Level 2 coarse-coarse lattice: " << clatt2 << std::endl; + + GridCartesian *CoarseCoarse4d = SpaceTimeGrid::makeFourDimGrid(clatt2, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); + GridCartesian *CoarseCoarse5d = SpaceTimeGrid::makeFiveDimGrid(1,CoarseCoarse4d); + + std::vector seeds4({1,2,3,4}); + std::vector seeds5({5,6,7,8}); + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + + LatticeGaugeField Umu(UGrid); + std::cout << GridLogMessage << "Reading gauge field" << std::endl; + FieldMetaData header; + std::string file("/ccs/home/poare/ckpoint_lat.1000"); + NerscIO::readConfiguration(Umu,header,file); + + RealD b_ = 1.5; + RealD c_ = 0.5; + 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; + typedef ShiftedPVdagMLinearOperator ShiftedPVdagM_t; + typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef LittleDiracOperator::CoarseVector CoarseVector; + typedef Aggregation Subspace; + typedef MGPreconditioner TwoLevelMG; + + PVdagM_t PVdagM(Ddwf,Dpv); + ShiftedPVdagM_t ShiftedPVdagM(0.01,Ddwf,Dpv); + + NextToNearestStencilGeometry5D geom(Coarse5d); + + // Subspace cache: save after generation, reload on subsequent runs to skip expensive setup. + // Set SUBSPACE_FILE to override the default path. + 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")); + + // Check if subspace file exists (boss rank checks, result broadcast via GlobalSum). + 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); + // Insurance: GLOBAL (whole-lattice) orthonormalise, matching CreateSubspaceGCR + // (Aggregates.h:196), in case the cached file predates it. Span-preserving + // and globally orthonormal -- NOT the block Orthogonalise() below, which would + // defeat the raw-null discipline (runMG promotes the RAW subspace to build L2; + // block-GS here -> psi_coarse = e_k). The raw copy in runMG happens AFTER this. + AggregatesGCR.GlobalOrthonormalise(); + // AggregatesGCR.Orthogonalise(); + std::cout << GridLogMessage << "Subspace loaded, globally orthonormalised (raw block basis preserved)." << std::endl; + } else { + std::cout << GridLogMessage << "*** GCR subspace generation ***" << std::endl; + AggregatesGCR.CreateSubspaceGCR(RNG5,PVdagM,nbasis); + std::cout << GridLogMessage << "Subspace generation: PVdagM operator uses:" << std::endl; + PVdagM.getApplications(); + PVdagM.clear(); + saveSubspace(AggregatesGCR.subspace, subspace_file); + std::cout << GridLogMessage << "Subspace saved to: " << subspace_file << std::endl; + } + + runMG( + FGrid, + Coarse5d, + CoarseCoarse5d, + geom, + PVdagM, + ShiftedPVdagM, + AggregatesGCR + ); + + std::cout << GridLogMessage << "Done" << std::endl; + Grid_finalize(); + return 0; +} diff --git a/examples/Example_pvdagm_3level_madj.cc b/examples/Example_pvdagm_3level_madj.cc new file mode 100644 index 000000000..53ad40b9e --- /dev/null +++ b/examples/Example_pvdagm_3level_madj.cc @@ -0,0 +1,625 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./examples/Example_pvdagm_3level.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 + +#include +#include +#include + +using namespace std; +using namespace Grid; + +template void readFile(T& out, std::string const fname){ + #ifdef HAVE_LIME + std::cout << Grid::GridLogMessage << "Reading: " << fname << std::endl; + Grid::emptyUserRecord record; + Grid::ScidacReader SR; + SR.open(fname); + SR.readScidacFieldRecord(out, record); + SR.close(); + #endif +} + +template void writeFile(T& in, std::string const fname){ + #ifdef HAVE_LIME + std::cout << Grid::GridLogMessage << "Writing: " << fname << std::endl; + Grid::emptyUserRecord record; + Grid::ScidacWriter SW(in.Grid()->IsBoss()); + SW.open(fname); + SW.writeScidacFieldRecord(in, record); + SW.close(); + #endif +} + +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; + int nApp; + int nAppDag; +public: + PVdagMLinearOperator(Matrix &Mat,Matrix &PV): _Mat(Mat),_PV(PV), nApp(0), nAppDag(0) {}; + + 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); + nApp++; + } + void AdjOp (const Field &in, Field &out){ + Field tmp(in.Grid()); + _PV.M(in,tmp); + _Mat.Mdag(tmp,out); + nAppDag++; + } + void clear() { nApp = 0; nAppDag = 0; } + void getApplications() { + std::cout << GridLogMessage << "# applications of PVdagM: " << nApp << std::endl; + std::cout << GridLogMessage << "# applications of PVdagM^dag: " << nAppDag << std::endl; + std::cout << GridLogMessage << "# applications total: " << nApp + nAppDag << std::endl; + } + 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 MdagPVLinearOperator : public LinearOperatorBase { + Matrix &_Mat; + Matrix &_PV; +public: + MdagPVLinearOperator(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()); + _PV.M(in,tmp); + _Mat.Mdag(tmp,out); + } + void AdjOp (const Field &in, Field &out){ + Field tmp(in.Grid()); + _Mat.M(in,tmp); + _PV.Mdag(tmp,out); + } + void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ + 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; + RealD shift; +public: + 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 MGPreconditioner : public LinearFunction< Lattice > { +public: + using LinearFunction >::operator(); + + typedef Aggregation Aggregates; + typedef typename Aggregation::FineField FineField; + typedef typename Aggregation::CoarseVector CoarseVector; + typedef typename Aggregation::CoarseMatrix CoarseMatrix; + typedef LinearOperatorBase FineOperator; + typedef LinearFunction FineSmoother; + typedef LinearOperatorBase CoarseOperator; + typedef LinearFunction CoarseSolver; + + Aggregates & _Aggregates; + FineOperator & _FineOperator; + FineSmoother & _PreSmoother; + FineSmoother & _PostSmoother; + CoarseOperator & _CoarseOperator; + CoarseSolver & _CoarseSolve; + CoarseSolver & _CoarseGuesser; + std::string _name; + + int level; void Level(int lv) {level = lv; }; + + MGPreconditioner(Aggregates &Agg, + FineOperator &Fine, + FineSmoother &PreSmoother, + FineSmoother &PostSmoother, + CoarseOperator &CoarseOperator_, + CoarseSolver &CoarseSolve_, + CoarseSolver &CoarseGuesser_, + std::string name = std::string("unnamed")) + : _Aggregates(Agg), + _FineOperator(Fine), + _PreSmoother(PreSmoother), + _PostSmoother(PostSmoother), + _CoarseOperator(CoarseOperator_), + _CoarseSolve(CoarseSolve_), + _CoarseGuesser(CoarseGuesser_), + _name(name), + level(1) { } + + virtual void operator()(const FineField &in, FineField & out) + { + GridBase *CoarseGrid = _Aggregates.CoarseGrid; + CoarseVector Csrc(CoarseGrid); + CoarseVector Csol(CoarseGrid); + FineField vec1(in.Grid()); + FineField vec2(in.Grid()); + + double t; + out = Zero(); + t=-usecond(); + _PreSmoother(in,out); + t+=usecond(); + std::cout< +class ShiftedLinearOperator : public LinearOperatorBase { + LinearOperatorBase &_Op; + RealD shift; +public: + ShiftedLinearOperator(RealD _shift, LinearOperatorBase &Op) : shift(_shift), _Op(Op) {} + 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); } +}; + +template +void runMG( + GridCartesian *FGrid, + GridCartesian *Coarse5d, + GridCartesian *CoarseCoarse5d, + NextToNearestStencilGeometry5D geom, + PVdagM_t &PVdagM, + ShiftedPVdagM_t &ShiftedPVdagM, + Subspace &AggregatesPD +) { + std::vector subspace = AggregatesPD.subspace; + assert((int)subspace.size() == NB); + const int nbasis = NB; + const int cb = 0; + + CoarseVector c_src(Coarse5d); + CoarseVector c_res(Coarse5d); + Complex one(1.0); + + LatticeFermionD f_src(FGrid); + LatticeFermionD f_res(FGrid); + + TrivialPrecon simpleC; + TrivialPrecon simple_fine; + + ////////////////////////////////////////////////////////////////////// + // Level 0→1: coarsen PVdagM, build LinOpCoarse + ////////////////////////////////////////////////////////////////////// + LittleDiracOperator LittleDiracOpPV(geom, FGrid, Coarse5d); + LittleDiracOpPV.CoarsenOperator(PVdagM, AggregatesPD); + + NonHermitianLinearOperator LinOpCoarse(LittleDiracOpPV); + + ////////////////////////////////////////////////////////////////////// + // Baseline: plain PGCR on LinOpCoarse (reference for comparison) + ////////////////////////////////////////////////////////////////////// + std::cout< L2PGCR_baseline(3.0e-2,1100,LinOpCoarse,simpleC,10,10); + L2PGCR_baseline.Level(2); + c_src = one; + c_res = Zero(); + L2PGCR_baseline(c_src,c_res); + PVdagM.getApplications(); + PVdagM.clear(); + + ////////////////////////////////////////////////////////////////////// + // psi_coarse: coarse projections of pre-GS fine null vectors. + // These are the Level 1 near-null vectors, promoted from Level 0. + // Used as the aggregation basis for Level 1→2 coarsening. + ////////////////////////////////////////////////////////////////////// + std::vector psi_coarse(nbasis, Coarse5d); + for (int k = 0; k < nbasis; k++) + AggregatesPD.ProjectToSubspace(psi_coarse[k], subspace[k]); + + ////////////////////////////////////////////////////////////////////// + // Diagnostics: W (fine projected matrix) and C (Galerkin check) + ////////////////////////////////////////////////////////////////////// + { + Eigen::MatrixXcd W = Eigen::MatrixXcd::Zero(nbasis, nbasis); + LatticeFermion ftmp(FGrid); + for (int j = 0; j < nbasis; j++) { + PVdagM.Op(subspace[j], ftmp); + for (int i = 0; i < nbasis; i++) + W(i,j) = TensorRemove(innerProduct(subspace[i], ftmp)); + } + RealD normW = W.norm(); + std::cout << GridLogMessage << "Fine projected matrix ||W|| = " << normW << std::endl; + + Eigen::MatrixXcd C = Eigen::MatrixXcd::Zero(nbasis, nbasis); + CoarseVector Ac(Coarse5d); + for (int l = 0; l < nbasis; l++) { + LinOpCoarse.Op(psi_coarse[l], Ac); + for (int k = 0; k < nbasis; k++) + C(k,l) = TensorRemove(innerProduct(psi_coarse[k], Ac)); + } + RealD normC = C.norm(); + RealD normCmCdag = (C - C.adjoint()).norm(); + std::cout << GridLogMessage << "Coarse null matrix ||C|| = " << normC << std::endl; + std::cout << GridLogMessage << "Coarse null matrix ||C - C†||/||C|| = " << normCmCdag/normC << std::endl; + std::cout << GridLogMessage << "Galerkin check ||C||/||W|| = " << normC/normW << std::endl; + } + + ////////////////////////////////////////////////////////////////////// + // Level 1→2: set up aggregation using psi_coarse as subspace. + // Block factor 2,2,3,2 (removes odd local sublattice in z given MPI + // geometry 3×6×4×4 where z-local at Level 1 is 6). + // psi_coarse are assigned directly; CoarsenOperator performs + // block-GS orthogonalisation before building LinOpCoarseCoarse. + ////////////////////////////////////////////////////////////////////// + // innerProduct(CoarseSiteObj, CoarseSiteObj) returns iScalar, so CComplex + // for the L1→L2 level must be iScalar, not vTComplex. + typedef typename CoarseVector::vector_object CoarseSiteObj; + typedef iScalar vTTComplex; + typedef GeneralCoarsenedMatrix LittleDiracOperatorL2; + typedef typename LittleDiracOperatorL2::CoarseVector CoarseCoarseVector; + typedef Aggregation SubspaceL2; + typedef MGPreconditioner L1to2MG; + + SubspaceL2 AggregatesL2(CoarseCoarse5d, Coarse5d, cb); + for (int k = 0; k < nbasis; k++) + AggregatesL2.subspace[k] = psi_coarse[k]; + + NextToNearestStencilGeometry5D geom2(CoarseCoarse5d); + LittleDiracOperatorL2 LittleDiracOpL2(geom2, Coarse5d, CoarseCoarse5d); + LittleDiracOpL2.CoarsenOperator(LinOpCoarse, AggregatesL2); + + NonHermitianLinearOperator LinOpCC(LittleDiracOpL2); + + ////////////////////////////////////////////////////////////////////// + // Level 2 solver: plain GCR, no further coarsening + ////////////////////////////////////////////////////////////////////// + TrivialPrecon simpleCC; + // L3PGCR is an inner solver inside the L1→2 V-cycle; does not need to converge + // to fine-grid precision. Loose tolerance (3e-2) and large restart (64) to allow + // the Krylov space to span enough of the near-null spectrum of LinOpCC per cycle. + PrecGeneralisedConjugateResidualNonHermitian L3PGCR(1.0e-4,5,LinOpCC,simpleCC,64,64); + L3PGCR.Level(3); + + ////////////////////////////////////////////////////////////////////// + // Coarse-level GCR smoother for Level 1→2 V-cycle. + // Mirrors fine-grid SmootherGCR: shifted operator + fixed step count. + // coarse_smoother_shift and coarse_smoother_nstep are the tuning knobs. + ////////////////////////////////////////////////////////////////////// + RealD coarse_smoother_shift = 0.0; + int coarse_smoother_nstep = 8; + if(getenv("coarse_smoother_shift")) coarse_smoother_shift = atof(getenv("coarse_smoother_shift")); + if(getenv("coarse_smoother_nstep")) coarse_smoother_nstep = atoi(getenv("coarse_smoother_nstep")); + + ShiftedLinearOperator ShiftedLinOpCoarse(coarse_smoother_shift, LinOpCoarse); + PrecGeneralisedConjugateResidualNonHermitian CoarseSmootherGCR(0.0,1,ShiftedLinOpCoarse,simpleC,coarse_smoother_nstep,coarse_smoother_nstep); + CoarseSmootherGCR.Level(2); + + ////////////////////////////////////////////////////////////////////// + // Level 1→2 V-cycle preconditioner. + ////////////////////////////////////////////////////////////////////// + L1to2MG L1to2Precon(AggregatesL2, + LinOpCoarse, + simpleC, // no pre-smoother (matches fine-grid setup) + CoarseSmootherGCR, // post-smoother: 12 GCR steps + LinOpCC, + L3PGCR, + simpleCC, + std::string("LinOpC")); + + ////////////////////////////////////////////////////////////////////// + // Standalone Level 1 two-level solve test. + // Compare against plain PGCR baseline above. + ////////////////////////////////////////////////////////////////////// + std::cout< L2MGsolver(3.0e-2,200,LinOpCoarse,L1to2Precon,16,16); + L2MGsolver.Level(2); + c_res = Zero(); + L2MGsolver(c_src,c_res); + + std::cout << GridLogMessage << "Level 1 two-level test: PVdagM operator uses:" << std::endl; + PVdagM.getApplications(); + PVdagM.clear(); + + ////////////////////////////////////////////////////////////////////// + // Full three-level outer solve + ////////////////////////////////////////////////////////////////////// + std::cout< SmootherGCR(0.00,1,ShiftedPVdagM,simple_fine,16,16); + SmootherGCR.Level(1); + + f_src = one; + + // Pre-smoother: none (TrivialPrecon); post-smoother: shifted PGCR. + // Coarse solver: L2MGsolver (PGCR preconditioned by Level 1→2 V-cycle). + TwoLevelMG ThreeLevelPrecon(AggregatesPD, + PVdagM, + simple_fine, + SmootherGCR, + LinOpCoarse, + L2MGsolver, + simpleC, + std::string("PVdagM")); + + PrecGeneralisedConjugateResidualNonHermitian L1PGCR(1.0e-8,1000,PVdagM,ThreeLevelPrecon,16,16); + L1PGCR.Level(1); + + f_res = Zero(); + L1PGCR(f_src,f_res); + + std::cout << GridLogMessage << "Three-level outer solve: PVdagM operator uses:" << std::endl; + PVdagM.getApplications(); + PVdagM.clear(); +} + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + const int Ls = 24; + RealD M5 = 1.8; + RealD b = 1.5; + RealD c = 0.5; + RealD mass = 0.00078; + if ( getenv("MASS") ) mass = atof(getenv("MASS")); + + const int nbasis = 60; + + std::cout << GridLogMessage << "Mass: " << mass << ", Ls: " << Ls << ", b=" << b << ", c=" << c << std::endl; + std::cout << GridLogMessage << "nbasis: " << nbasis << std::endl; + + 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); + + // Level 1 coarse grid: block 2^4 from fine (48×48×48×96 → 24×24×24×48, Ls=1) + Coordinate clatt = lat_size; + for (int d = 0; d < 4; d++) clatt[d] /= 2; + std::cout << GridLogMessage << "Level 1 coarse lattice: " << clatt << std::endl; + + GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); + GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d); + + // Level 2 coarse-coarse grid: block 2,2,3,3 from Level 1 (24×24×24×48 → 12×12×8×16, Ls=1). + // MPI geometry 3.6.4.4 (288 ranks): fine local {16,8,12,24}. + // Level 1 local {8,4,6,12}; Level 2 local {4,2,2,4}. + // z blocked by 3: z-Level1-local=6; 6/3=2 (even), 6/2=3 (odd) → must use 3. + // t blocked by 3: t-Level1-local=12; 12/3=4 divisible by Nsimd=4 (gen-simd-width=64). + // t-block=2 gives t2-local=6, 6 mod 4 ≠ 0, fails Grid SIMD assertion. ✓ + // With {4,2,2,4}: Nsimd=4 goes into x or t (both =4). ✓ + Coordinate clatt2 = clatt; + clatt2[0] /= 2; + clatt2[1] /= 2; + clatt2[2] /= 3; + clatt2[3] /= 3; + std::cout << GridLogMessage << "Level 2 coarse-coarse lattice: " << clatt2 << std::endl; + + GridCartesian *CoarseCoarse4d = SpaceTimeGrid::makeFourDimGrid(clatt2, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); + GridCartesian *CoarseCoarse5d = SpaceTimeGrid::makeFiveDimGrid(1,CoarseCoarse4d); + + std::vector seeds4({1,2,3,4}); + std::vector seeds5({5,6,7,8}); + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + + LatticeGaugeField Umu(UGrid); + std::cout << GridLogMessage << "Reading gauge field" << std::endl; + FieldMetaData header; + std::string file("/ccs/home/poare/ckpoint_lat.1000"); + NerscIO::readConfiguration(Umu,header,file); + + RealD b_ = 1.5; + RealD c_ = 0.5; + + RealD madj = 1.0; + if ( getenv("MADJ") ) madj=atof(getenv("MADJ")); + std::cout << "PV mass set to "< PVdagM_t; + typedef ShiftedPVdagMLinearOperator ShiftedPVdagM_t; + typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef LittleDiracOperator::CoarseVector CoarseVector; + typedef Aggregation Subspace; + typedef MGPreconditioner TwoLevelMG; + + PVdagM_t PVdagM(Ddwf,Dpv); + ShiftedPVdagM_t ShiftedPVdagM(0.00,Ddwf,Dpv); + + NextToNearestStencilGeometry5D geom(Coarse5d); + + // Subspace cache: save after generation, reload on subsequent runs to skip expensive setup. + // Set SUBSPACE_FILE to override the default path. + 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")); + + // Check if subspace file exists (boss rank checks, result broadcast via GlobalSum). + 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); + // Re-orthogonalise after loading to ensure block-GS condition holds. + // AggregatesGCR.Orthogonalise(); + std::cout << GridLogMessage << "Subspace loaded and re-orthogonalised." << std::endl; + } else { + std::cout << GridLogMessage << "*** GCR subspace generation ***" << std::endl; + AggregatesGCR.CreateSubspaceGCR(RNG5,PVdagM,nbasis); + std::cout << GridLogMessage << "Subspace generation: PVdagM operator uses:" << std::endl; + PVdagM.getApplications(); + PVdagM.clear(); + saveSubspace(AggregatesGCR.subspace, subspace_file); + std::cout << GridLogMessage << "Subspace saved to: " << subspace_file << std::endl; + } + + runMG( + FGrid, + Coarse5d, + CoarseCoarse5d, + geom, + PVdagM, + ShiftedPVdagM, + AggregatesGCR + ); + + std::cout << GridLogMessage << "Done" << std::endl; + Grid_finalize(); + return 0; +} diff --git a/examples/Example_pvdagm_4level.cc b/examples/Example_pvdagm_4level.cc new file mode 100644 index 000000000..2d8d67964 --- /dev/null +++ b/examples/Example_pvdagm_4level.cc @@ -0,0 +1,798 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./examples/Example_pvdagm_3level.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 + +#include +#include +#include + +using namespace std; +using namespace Grid; + +template void readFile(T& out, std::string const fname){ + #ifdef HAVE_LIME + std::cout << Grid::GridLogMessage << "Reading: " << fname << std::endl; + Grid::emptyUserRecord record; + Grid::ScidacReader SR; + SR.open(fname); + SR.readScidacFieldRecord(out, record); + SR.close(); + #endif +} + +template void writeFile(T& in, std::string const fname){ + #ifdef HAVE_LIME + std::cout << Grid::GridLogMessage << "Writing: " << fname << std::endl; + Grid::emptyUserRecord record; + Grid::ScidacWriter SW(in.Grid()->IsBoss()); + SW.open(fname); + SW.writeScidacFieldRecord(in, record); + SW.close(); + #endif +} + +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; + int nApp; + int nAppDag; +public: + PVdagMLinearOperator(Matrix &Mat,Matrix &PV): _Mat(Mat),_PV(PV), nApp(0), nAppDag(0) {}; + + 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); + nApp++; + } + void AdjOp (const Field &in, Field &out){ + Field tmp(in.Grid()); + _PV.M(in,tmp); + _Mat.Mdag(tmp,out); + nAppDag++; + } + void clear() { nApp = 0; nAppDag = 0; } + void getApplications() { + std::cout << GridLogMessage << "# applications of PVdagM: " << nApp << std::endl; + std::cout << GridLogMessage << "# applications of PVdagM^dag: " << nAppDag << std::endl; + std::cout << GridLogMessage << "# applications total: " << nApp + nAppDag << std::endl; + } + 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 MdagPVLinearOperator : public LinearOperatorBase { + Matrix &_Mat; + Matrix &_PV; +public: + MdagPVLinearOperator(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()); + _PV.M(in,tmp); + _Mat.Mdag(tmp,out); + } + void AdjOp (const Field &in, Field &out){ + Field tmp(in.Grid()); + _Mat.M(in,tmp); + _PV.Mdag(tmp,out); + } + void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ + 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; + RealD shift; +public: + 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); + } +}; + +// Lüscher deflated guesser (arXiv:0706.2298 Sec A.3) for a non-Hermitian solve. +// C_{st} = ; guess = sum_s c_s psi[s] where c = C^{-1} psi† src. +template +class LuscherGuesser : public LinearFunction { + const std::vector ψ + Eigen::MatrixXcd C_inv; +public: + using LinearFunction::operator(); + LuscherGuesser(const std::vector &psi_, const Eigen::MatrixXcd &Cinv_) + : psi(psi_), C_inv(Cinv_) {} + virtual void operator()(const Field &src, Field &guess) { + int N = psi.size(); + Eigen::VectorXcd b(N); + for (int t = 0; t < N; t++) + b(t) = TensorRemove(innerProduct(psi[t], src)); + Eigen::VectorXcd c = C_inv * b; + guess = Zero(); + for (int s = 0; s < N; s++) + guess += ComplexD(c(s)) * psi[s]; + } +}; + +template +class MGPreconditioner : public LinearFunction< Lattice > { +public: + using LinearFunction >::operator(); + + typedef Aggregation Aggregates; + typedef typename Aggregation::FineField FineField; + typedef typename Aggregation::CoarseVector CoarseVector; + typedef typename Aggregation::CoarseMatrix CoarseMatrix; + typedef LinearOperatorBase FineOperator; + typedef LinearFunction FineSmoother; + typedef LinearOperatorBase CoarseOperator; + typedef LinearFunction CoarseSolver; + + Aggregates & _Aggregates; + FineOperator & _FineOperator; + FineSmoother & _PreSmoother; + FineSmoother & _PostSmoother; + CoarseOperator & _CoarseOperator; + CoarseSolver & _CoarseSolve; + CoarseSolver & _CoarseGuesser; + + int level; void Level(int lv) {level = lv; }; + + MGPreconditioner(Aggregates &Agg, + FineOperator &Fine, + FineSmoother &PreSmoother, + FineSmoother &PostSmoother, + CoarseOperator &CoarseOperator_, + CoarseSolver &CoarseSolve_, + CoarseSolver &CoarseGuesser_) + : _Aggregates(Agg), + _FineOperator(Fine), + _PreSmoother(PreSmoother), + _PostSmoother(PostSmoother), + _CoarseOperator(CoarseOperator_), + _CoarseSolve(CoarseSolve_), + _CoarseGuesser(CoarseGuesser_), + level(1) { } + + virtual void operator()(const FineField &in, FineField & out) + { + GridBase *CoarseGrid = _Aggregates.CoarseGrid; + CoarseVector Csrc(CoarseGrid); + CoarseVector Csol(CoarseGrid); + FineField vec1(in.Grid()); + FineField vec2(in.Grid()); + + double t; + out = Zero(); + t=-usecond(); + _PreSmoother(in,out); + t+=usecond(); + std::cout< +class ShiftedLinearOperator : public LinearOperatorBase { + LinearOperatorBase &_Op; + RealD shift; +public: + ShiftedLinearOperator(RealD _shift, LinearOperatorBase &Op) : shift(_shift), _Op(Op) {} + 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); } +}; + +template +void runMG( + GridCartesian *FGrid, + GridCartesian *Coarse5d, + GridCartesian *CoarseCoarse5d, + GridCartesian *CoarseCoarseCoarse5d, + NextToNearestStencilGeometry5D geom, + PVdagM_t &PVdagM, + ShiftedPVdagM_t &ShiftedPVdagM, + Subspace &AggregatesPD +) { + std::vector subspace = AggregatesPD.subspace; + assert((int)subspace.size() == NB); + const int nbasis = NB; + const int cb = 0; + + CoarseVector c_src(Coarse5d); + CoarseVector c_res(Coarse5d); + Complex one(1.0); + + LatticeFermionD f_src(FGrid); + LatticeFermionD f_res(FGrid); + + TrivialPrecon simpleC; + TrivialPrecon simple_fine; + + ////////////////////////////////////////////////////////////////////// + // Level 0→1: coarsen PVdagM, build LinOpCoarse + ////////////////////////////////////////////////////////////////////// + LittleDiracOperator LittleDiracOpPV(geom, FGrid, Coarse5d); + LittleDiracOpPV.CoarsenOperator(PVdagM, AggregatesPD); + + NonHermitianLinearOperator LinOpCoarse(LittleDiracOpPV); + + ////////////////////////////////////////////////////////////////////// + // Baseline: plain PGCR on LinOpCoarse (reference for comparison) + ////////////////////////////////////////////////////////////////////// + std::cout< L2PGCR_baseline(3.0e-2,1100,LinOpCoarse,simpleC,10,10); + L2PGCR_baseline.Level(2); + L2PGCR_baseline.Name("Cbaseline"); + c_src = one; + c_res = Zero(); + L2PGCR_baseline(c_src,c_res); + + ////////////////////////////////////////////////////////////////////// + // psi_coarse: coarse projections of pre-GS fine null vectors. + // These are the Level 1 near-null vectors, promoted from Level 0. + // Used as the aggregation basis for Level 1→2 coarsening. + ////////////////////////////////////////////////////////////////////// + std::vector psi_coarse(nbasis, Coarse5d); + for (int k = 0; k < nbasis; k++) + AggregatesPD.ProjectToSubspace(psi_coarse[k], subspace[k]); + + ////////////////////////////////////////////////////////////////////// + // Diagnostics: W (fine projected matrix) and C (Galerkin check) + ////////////////////////////////////////////////////////////////////// + { + Eigen::MatrixXcd W = Eigen::MatrixXcd::Zero(nbasis, nbasis); + LatticeFermion ftmp(FGrid); + for (int j = 0; j < nbasis; j++) { + PVdagM.Op(subspace[j], ftmp); + for (int i = 0; i < nbasis; i++) + W(i,j) = TensorRemove(innerProduct(subspace[i], ftmp)); + } + RealD normW = W.norm(); + std::cout << GridLogMessage << "Fine projected matrix ||W|| = " << normW << std::endl; + + Eigen::MatrixXcd C = Eigen::MatrixXcd::Zero(nbasis, nbasis); + CoarseVector Ac(Coarse5d); + for (int l = 0; l < nbasis; l++) { + LinOpCoarse.Op(psi_coarse[l], Ac); + for (int k = 0; k < nbasis; k++) + C(k,l) = TensorRemove(innerProduct(psi_coarse[k], Ac)); + } + RealD normC = C.norm(); + RealD normCmCdag = (C - C.adjoint()).norm(); + std::cout << GridLogMessage << "Coarse null matrix ||C|| = " << normC << std::endl; + std::cout << GridLogMessage << "Coarse null matrix ||C - C†||/||C|| = " << normCmCdag/normC << std::endl; + std::cout << GridLogMessage << "Galerkin check ||C||/||W|| = " << normC/normW << std::endl; + } + + ////////////////////////////////////////////////////////////////////// + // Level 1→2: set up aggregation using psi_coarse as subspace. + // Block factor 2,2,3,2 (removes odd local sublattice in z given MPI + // geometry 3×6×4×4 where z-local at Level 1 is 6). + // psi_coarse are assigned directly; CoarsenOperator performs + // block-GS orthogonalisation before building LinOpCoarseCoarse. + ////////////////////////////////////////////////////////////////////// + // innerProduct(CoarseSiteObj, CoarseSiteObj) returns iScalar, so CComplex + // for the L1→L2 level must be iScalar, not vTComplex. + typedef typename CoarseVector::vector_object CoarseSiteObj; + typedef iScalar vTTComplex; + typedef GeneralCoarsenedMatrix LittleDiracOperatorL2; + typedef typename LittleDiracOperatorL2::CoarseVector CoarseCoarseVector; + typedef Aggregation SubspaceL2; + typedef MGPreconditioner L1to2MG; + + SubspaceL2 AggregatesL2(CoarseCoarse5d, Coarse5d, cb); + for (int k = 0; k < nbasis; k++) + AggregatesL2.subspace[k] = psi_coarse[k]; + + NextToNearestStencilGeometry5D geom2(CoarseCoarse5d); + LittleDiracOperatorL2 LittleDiracOpL2(geom2, Coarse5d, CoarseCoarse5d); + LittleDiracOpL2.CoarsenOperator(LinOpCoarse, AggregatesL2); + + NonHermitianLinearOperator LinOpCC(LittleDiracOpL2); + + TrivialPrecon simpleCC; + + ////////////////////////////////////////////////////////////////////// + // Lüscher deflation guesser for L3PGCR. + // Step 1: project psi_coarse[k] (promoted fine null vectors) to + // CoarseCoarseVector space — these cover the zero-momentum + // component of the near-null space of LinOpCC. + // Step 2: breed Nextra additional null vectors directly on LinOpCC + // using GCR with random sources — these pick up near-null + // modes at all spatial frequencies not spanned by step 1. + // Step 3: build C_{st} = over the + // full augmented basis and invert directly via Eigen LU. + ////////////////////////////////////////////////////////////////////// + std::vector psi_cc(nbasis, CoarseCoarse5d); + for (int k = 0; k < nbasis; k++) + AggregatesL2.ProjectToSubspace(psi_cc[k], psi_coarse[k]); + + { + int Nextra = nbasis; // breed as many extra as we have promoted ones + if ( getenv("CC_NEXTRA") ) Nextra = atoi(getenv("CC_NEXTRA")); + GridParallelRNG RNG_CC(CoarseCoarse5d); + RNG_CC.SeedFixedIntegers({11,13,17,19}); + PrecGeneralisedConjugateResidualNonHermitian + nullGCR(1e-2, 200, LinOpCC, simpleCC, 32, 32); + CoarseCoarseVector tmp(CoarseCoarse5d); + for (int k = 0; k < Nextra; k++) { + CoarseCoarseVector src(CoarseCoarse5d); + gaussian(RNG_CC, src); + tmp = Zero(); + nullGCR(src, tmp); + psi_cc.push_back(tmp); + } + std::cout << GridLogMessage << "LinOpCC deflation basis: " << nbasis + << " promoted + " << Nextra << " bred = " << psi_cc.size() << " total" << std::endl; + } + + const int Naug = psi_cc.size(); + Eigen::MatrixXcd Ccc = Eigen::MatrixXcd::Zero(Naug, Naug); + { + CoarseCoarseVector Acc(CoarseCoarse5d); + for (int l = 0; l < Naug; l++) { + LinOpCC.Op(psi_cc[l], Acc); + for (int k = 0; k < Naug; k++) + Ccc(k,l) = TensorRemove(innerProduct(psi_cc[k], Acc)); + } + } + { + RealD normCcc = Ccc.norm(); + RealD normCccmCdag = (Ccc - Ccc.adjoint()).norm(); + std::cout << GridLogMessage << "Coarse-coarse deflation matrix ||Ccc|| = " << normCcc << std::endl; + std::cout << GridLogMessage << "Coarse-coarse deflation matrix ||Ccc-Ccc†||/||Ccc|| = " << normCccmCdag/normCcc << std::endl; + } + Eigen::MatrixXcd Ccc_inv = Ccc.inverse(); + LuscherGuesser CCDeflGuesser(psi_cc, Ccc_inv); + + ////////////////////////////////////////////////////////////////////// + // Level 2→3: coarsen LinOpCC using the RAW promoted psi_cc as aggregation + // to build the Level 4 (coarse-coarse-coarse) operator. + // psi_cc[0..nbasis-1] are the coarse-coarse near-null vectors, projected + // from the RAW psi_coarse (themselves projected from the RAW fine null + // vectors) -- the pre-block-GS chain the whole construction depends on. + // CoarsenOperator block-GS orthogonalises AggregatesL3.subspace IN PLACE, + // so assign COPIES of psi_cc and keep psi_cc itself raw. + // + // Tensor depth deepens once more: innerProduct(CoarseCoarseSiteObj,...) returns + // iScalar, so CComplex for the L2→L3 level is iScalar>. + ////////////////////////////////////////////////////////////////////// + typedef typename CoarseCoarseVector::vector_object CoarseCoarseSiteObj; + typedef iScalar vTTTComplex; + typedef GeneralCoarsenedMatrix LittleDiracOperatorL3; + typedef typename LittleDiracOperatorL3::CoarseVector CoarseCoarseCoarseVector; + typedef Aggregation SubspaceL3; + typedef MGPreconditioner L2to3MG; + + SubspaceL3 AggregatesL3(CoarseCoarseCoarse5d, CoarseCoarse5d, cb); + for (int k = 0; k < nbasis; k++) + AggregatesL3.subspace[k] = psi_cc[k]; // raw promoted; COPY, keeps psi_cc raw + + NextToNearestStencilGeometry5D geom3(CoarseCoarseCoarse5d); + LittleDiracOperatorL3 LittleDiracOpL3(geom3, CoarseCoarse5d, CoarseCoarseCoarse5d); + LittleDiracOpL3.CoarsenOperator(LinOpCC, AggregatesL3); // block-GS's AggregatesL3.subspace in place + + NonHermitianLinearOperator LinOpCCC(LittleDiracOpL3); + TrivialPrecon simpleCCC; + + ////////////////////////////////////////////////////////////////////// + // Level 4 bottom solve: GCR on a SHIFTED LinOpCCC. This is the one level + // with no IRS shift, and it is the most non-normal (coarsest) operator, so + // the bare bottom GCR wanders in a field of values that wraps the origin and + // its iteration count blows out (observed 5..54 iters to hit 0.2). Solving + // (A_ccc + l4_shift) instead slides the FoV off the origin; the correction is + // only ever a loose 0.2 approximation anyway, so the detuning is free. + // l4_shift defaults to 0.0 => bare LinOpCCC, baseline unchanged until opted in. + ////////////////////////////////////////////////////////////////////// + RealD l4_shift = 0.0; + if(getenv("l4_shift")) l4_shift = atof(getenv("l4_shift")); + std::cout << GridLogMessage << "PARAM l4_shift = " << l4_shift << std::endl; + + ShiftedLinearOperator ShiftedLinOpCCC(l4_shift, LinOpCCC); + PrecGeneralisedConjugateResidualNonHermitian L4PGCR(1.0e-1,200,ShiftedLinOpCCC,simpleCCC,16,16); + L4PGCR.Level(4); + L4PGCR.Name("CCCouter"); + + ////////////////////////////////////////////////////////////////////// + // Level 2→3 V-cycle: depth-2 SHIFTED smoother on LinOpCC + L4 bottom. + // The shift slides the coarse-coarse field of values off the origin so a + // 2-step smoother has something to bite on a non-normal operator (IRS idea). + ////////////////////////////////////////////////////////////////////// + RealD cc_smoother_shift = 0.01; + int cc_smoother_nstep = 2; + if(getenv("cc_smoother_shift")) cc_smoother_shift = atof(getenv("cc_smoother_shift")); + if(getenv("cc_smoother_nstep")) cc_smoother_nstep = atoi(getenv("cc_smoother_nstep")); + + ShiftedLinearOperator ShiftedLinOpCC(cc_smoother_shift, LinOpCC); + PrecGeneralisedConjugateResidualNonHermitian + CoarseCoarseSmootherGCR(0.01,1,ShiftedLinOpCC,simpleCC,cc_smoother_nstep,cc_smoother_nstep); + CoarseCoarseSmootherGCR.SetZeroGuess(1); // smoother slot: caller zeroes guess + CoarseCoarseSmootherGCR.Level(3); + CoarseCoarseSmootherGCR.Name("CCsmoother"); + + L2to3MG L2to3Precon(AggregatesL3, + LinOpCC, + simpleCC, // no pre-smoother + CoarseCoarseSmootherGCR, // post-smoother: depth-2 shifted GCR + LinOpCCC, + L4PGCR, + simpleCCC); // trivial guesser at the bottom + + ////////////////////////////////////////////////////////////////////// + // Level 3 (coarse-coarse) solve: GCR preconditioned by the L2→L3 V-cycle. + // Replaces the plain L3PGCR of the 3-level build -- the coarse-coarse level + // is now smoothed shallowly and recursed rather than solved deeply. + ////////////////////////////////////////////////////////////////////// + PrecGeneralisedConjugateResidualNonHermitian L3MGsolver(1.0e-1,200,LinOpCC,L2to3Precon,16,16); + L3MGsolver.Level(3); + L3MGsolver.Name("CCouter"); + + ////////////////////////////////////////////////////////////////////// + // Coarse-level GCR smoother for Level 1→2 V-cycle. + // Mirrors fine-grid SmootherGCR: shifted operator + fixed step count. + // coarse_smoother_shift and coarse_smoother_nstep are the tuning knobs. + ////////////////////////////////////////////////////////////////////// + RealD coarse_smoother_shift = 0.01; + int coarse_smoother_nstep = 2; // depth-2 smoother on the coarse level + if(getenv("coarse_smoother_shift")) coarse_smoother_shift = atof(getenv("coarse_smoother_shift")); + if(getenv("coarse_smoother_nstep")) coarse_smoother_nstep = atoi(getenv("coarse_smoother_nstep")); + + ShiftedLinearOperator ShiftedLinOpCoarse(coarse_smoother_shift, LinOpCoarse); + PrecGeneralisedConjugateResidualNonHermitian CoarseSmootherGCR(0.01,1,ShiftedLinOpCoarse,simpleC,coarse_smoother_nstep,coarse_smoother_nstep); + CoarseSmootherGCR.SetZeroGuess(1); // smoother slot: caller zeroes guess + CoarseSmootherGCR.Level(2); + CoarseSmootherGCR.Name("Csmoother"); + + ////////////////////////////////////////////////////////////////////// + // Level 1→2 V-cycle preconditioner. + ////////////////////////////////////////////////////////////////////// + L1to2MG L1to2Precon(AggregatesL2, + LinOpCoarse, + simpleC, // no pre-smoother (matches fine-grid setup) + CoarseSmootherGCR, // post-smoother: depth-2 shifted GCR + LinOpCC, + L3MGsolver, // coarse-coarse solve is now the L2→L3 V-cycle + CCDeflGuesser); // Lüscher guesser: psi_cc C^{-1} psi_cc† + + ////////////////////////////////////////////////////////////////////// + // Standalone Level 1 two-level solve test. + // Compare against plain PGCR baseline above. + ////////////////////////////////////////////////////////////////////// + std::cout< L2MGsolver(3.0e-2,200,LinOpCoarse,L1to2Precon,16,16); + L2MGsolver.Level(2); + L2MGsolver.Name("Couter"); + c_res = Zero(); + L2MGsolver(c_src,c_res); + + std::cout << GridLogMessage << "Level 1 two-level test: PVdagM operator uses:" << std::endl; + PVdagM.getApplications(); + PVdagM.clear(); + + ////////////////////////////////////////////////////////////////////// + // Full three-level outer solve + ////////////////////////////////////////////////////////////////////// + std::cout< SmootherGCR(0.01,1,ShiftedPVdagM,simple_fine,16,16); + SmootherGCR.SetZeroGuess(1); // pre/post smoother slots zero their guess + SmootherGCR.Level(1); + SmootherGCR.Name("Fsmoother"); + + f_src = one; + + // Pre-smoother: none (TrivialPrecon); post-smoother: shifted PGCR. + // Coarse solver: L2MGsolver (PGCR preconditioned by Level 1→2 V-cycle). + TwoLevelMG ThreeLevelPrecon(AggregatesPD, + PVdagM, + simple_fine, + SmootherGCR, + LinOpCoarse, + L2MGsolver, + simpleC); + + PrecGeneralisedConjugateResidualNonHermitian L1PGCR(1.0e-8,1000,PVdagM,ThreeLevelPrecon,16,16); + L1PGCR.Level(1); + L1PGCR.Name("Fouter"); + + f_res = Zero(); + L1PGCR(f_src,f_res); + + std::cout << GridLogMessage << "Three-level outer solve: PVdagM operator uses:" << std::endl; + PVdagM.getApplications(); + PVdagM.clear(); +} + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + const int Ls = 24; + RealD M5 = 1.8; + RealD b = 1.5; + RealD c = 0.5; + RealD mass = 0.00078; + if ( getenv("MASS") ) mass = atof(getenv("MASS")); + + const int nbasis = 60; + + std::cout << GridLogMessage << "Mass: " << mass << ", Ls: " << Ls << ", b=" << b << ", c=" << c << std::endl; + std::cout << GridLogMessage << "nbasis: " << nbasis << std::endl; + + 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); + + // Level 1 coarse grid: block 2^4 from fine (48×48×48×96 → 24×24×24×48, Ls=1) + Coordinate clatt = lat_size; + for (int d = 0; d < 4; d++) clatt[d] /= 2; + std::cout << GridLogMessage << "Level 1 coarse lattice: " << clatt << std::endl; + + GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); + GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d); + + // Level 2 coarse-coarse grid: block 2,2,3,3 from Level 1 (24×24×24×48 → 12×12×8×16, Ls=1). + // MPI geometry 3.6.4.4 (288 ranks): fine local {16,8,12,24}. + // Level 1 local {8,4,6,12}; Level 2 local {4,2,2,4}. + // z blocked by 3: z-Level1-local=6; 6/3=2 (even), 6/2=3 (odd) → must use 3. + // t blocked by 3: t-Level1-local=12; 12/3=4 divisible by Nsimd=4 (gen-simd-width=64). + // t-block=2 gives t2-local=6, 6 mod 4 ≠ 0, fails Grid SIMD assertion. ✓ + // With {4,2,2,4}: Nsimd=4 goes into x or t (both =4). ✓ + Coordinate clatt2 = clatt; + clatt2[0] /= 2; + clatt2[1] /= 2; + clatt2[2] /= 3; + clatt2[3] /= 3; + std::cout << GridLogMessage << "Level 2 coarse-coarse lattice: " << clatt2 << std::endl; + + GridCartesian *CoarseCoarse4d = SpaceTimeGrid::makeFourDimGrid(clatt2, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); + GridCartesian *CoarseCoarse5d = SpaceTimeGrid::makeFiveDimGrid(1,CoarseCoarse4d); + + // Level 3 coarse-coarse-coarse grid: block clatt2 = {12,12,8,16} -> {6,12,8,8}. + // GEOMETRY (mpi 3.6.4.4, Nsimd=4 => SIMD layout {1,1,2,2}, factor 2 on z and t): + // every grid needs z-local and t-local EVEN. clatt2-local is {4,2,2,4}, so + // z-local=2 is already at its minimum even value and CANNOT be blocked (2->1 + // is odd and trips the SIMD assertion); y-local=2 would go to 1 (degenerate). + // Only x and t have room, so block {2,1,1,2}: clatt3 {6,12,8,8}, L4-local + // {2,2,2,2} -- all dims even and >=2. z stays unblocked by construction. + Coordinate clatt3 = clatt2; + clatt3[0] /= 2; // x: 12 -> 6 (x-local 4 -> 2) + // clatt3[1] (y) unblocked: y-local 2 -> blocking gives 1 (degenerate) + // clatt3[2] (z) unblocked: z-local 2 is SIMD-pinned even, cannot halve + clatt3[3] /= 2; // t: 16 -> 8 (t-local 4 -> 2) + std::cout << GridLogMessage << "Level 3 coarse-coarse-coarse lattice: " << clatt3 << std::endl; + + GridCartesian *CoarseCoarseCoarse4d = SpaceTimeGrid::makeFourDimGrid(clatt3, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); + GridCartesian *CoarseCoarseCoarse5d = SpaceTimeGrid::makeFiveDimGrid(1,CoarseCoarseCoarse4d); + + std::vector seeds4({1,2,3,4}); + std::vector seeds5({5,6,7,8}); + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + + LatticeGaugeField Umu(UGrid); + std::cout << GridLogMessage << "Reading gauge field" << std::endl; + FieldMetaData header; + std::string file("/ccs/home/poare/ckpoint_lat.1000"); + NerscIO::readConfiguration(Umu,header,file); + + RealD b_ = 1.5; + RealD c_ = 0.5; + 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; + typedef ShiftedPVdagMLinearOperator ShiftedPVdagM_t; + typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef LittleDiracOperator::CoarseVector CoarseVector; + typedef Aggregation Subspace; + typedef MGPreconditioner TwoLevelMG; + + PVdagM_t PVdagM(Ddwf,Dpv); + ShiftedPVdagM_t ShiftedPVdagM(0.01,Ddwf,Dpv); + + NextToNearestStencilGeometry5D geom(Coarse5d); + + // Subspace cache: save after generation, reload on subsequent runs to skip expensive setup. + // Set SUBSPACE_FILE to override the default path. + 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")); + + // Check if subspace file exists (boss rank checks, result broadcast via GlobalSum). + 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); + // Insurance: GLOBAL (whole-lattice) orthonormalise, in case the cached file + // predates the GlobalOrthonormalise() that CreateSubspaceGCR now applies + // (Aggregates.h:196). It is span-preserving and makes the vectors globally + // orthonormal -- it is NOT the block Orthogonalise() below, so it does NOT + // cause the psi_coarse->e_k trap. The RAW subspace copy in runMG happens + // AFTER this call, so the raw-null (pre-block-GS) discipline is preserved. + AggregatesGCR.GlobalOrthonormalise(); + // DO NOT block-orthogonalise here: runMG copies subspace[] as the RAW + // (pre-block-GS) basis and CoarsenOperator block-GS's it in place later. + // Orthogonalising now defeats the raw-null discipline (psi_coarse -> e_k) + // and poisons L2/L3/L4. See project_block_orthogonalise_leak. + // AggregatesGCR.Orthogonalise(); + std::cout << GridLogMessage << "Subspace loaded, globally orthonormalised (raw block basis preserved)." << std::endl; + } else { + std::cout << GridLogMessage << "*** GCR subspace generation ***" << std::endl; + AggregatesGCR.CreateSubspaceGCR(RNG5,PVdagM,nbasis); + std::cout << GridLogMessage << "Subspace generation: PVdagM operator uses:" << std::endl; + PVdagM.getApplications(); + PVdagM.clear(); + saveSubspace(AggregatesGCR.subspace, subspace_file); + std::cout << GridLogMessage << "Subspace saved to: " << subspace_file << std::endl; + } + + runMG( + FGrid, + Coarse5d, + CoarseCoarse5d, + CoarseCoarseCoarse5d, + geom, + PVdagM, + ShiftedPVdagM, + AggregatesGCR + ); + + std::cout << GridLogMessage << "Done" << std::endl; + Grid_finalize(); + return 0; +} diff --git a/examples/Example_pvdagm_5level.cc b/examples/Example_pvdagm_5level.cc new file mode 100644 index 000000000..115421cf1 --- /dev/null +++ b/examples/Example_pvdagm_5level.cc @@ -0,0 +1,954 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./examples/Example_pvdagm_5level.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 + +#include +#include +#include + +using namespace std; +using namespace Grid; + +template void readFile(T& out, std::string const fname){ + #ifdef HAVE_LIME + std::cout << Grid::GridLogMessage << "Reading: " << fname << std::endl; + Grid::emptyUserRecord record; + Grid::ScidacReader SR; + SR.open(fname); + SR.readScidacFieldRecord(out, record); + SR.close(); + #endif +} + +template void writeFile(T& in, std::string const fname){ + #ifdef HAVE_LIME + std::cout << Grid::GridLogMessage << "Writing: " << fname << std::endl; + Grid::emptyUserRecord record; + Grid::ScidacWriter SW(in.Grid()->IsBoss()); + SW.open(fname); + SW.writeScidacFieldRecord(in, record); + SW.close(); + #endif +} + +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; + int nApp; + int nAppDag; +public: + PVdagMLinearOperator(Matrix &Mat,Matrix &PV): _Mat(Mat),_PV(PV), nApp(0), nAppDag(0) {}; + + 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); + nApp++; + } + void AdjOp (const Field &in, Field &out){ + Field tmp(in.Grid()); + _PV.M(in,tmp); + _Mat.Mdag(tmp,out); + nAppDag++; + } + void clear() { nApp = 0; nAppDag = 0; } + void getApplications() { + std::cout << GridLogMessage << "# applications of PVdagM: " << nApp << std::endl; + std::cout << GridLogMessage << "# applications of PVdagM^dag: " << nAppDag << std::endl; + std::cout << GridLogMessage << "# applications total: " << nApp + nAppDag << std::endl; + } + 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 MdagPVLinearOperator : public LinearOperatorBase { + Matrix &_Mat; + Matrix &_PV; +public: + MdagPVLinearOperator(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()); + _PV.M(in,tmp); + _Mat.Mdag(tmp,out); + } + void AdjOp (const Field &in, Field &out){ + Field tmp(in.Grid()); + _Mat.M(in,tmp); + _PV.Mdag(tmp,out); + } + void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ + 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; + RealD shift; +public: + 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); + } +}; + +// Lüscher deflated guesser (arXiv:0706.2298 Sec A.3) for a non-Hermitian solve. +// C_{st} = ; guess = sum_s c_s psi[s] where c = C^{-1} psi† src. +template +class LuscherGuesser : public LinearFunction { + const std::vector ψ + Eigen::MatrixXcd C_inv; +public: + using LinearFunction::operator(); + LuscherGuesser(const std::vector &psi_, const Eigen::MatrixXcd &Cinv_) + : psi(psi_), C_inv(Cinv_) {} + virtual void operator()(const Field &src, Field &guess) { + int N = psi.size(); + Eigen::VectorXcd b(N); + for (int t = 0; t < N; t++) + b(t) = TensorRemove(innerProduct(psi[t], src)); + Eigen::VectorXcd c = C_inv * b; + guess = Zero(); + for (int s = 0; s < N; s++) + guess += ComplexD(c(s)) * psi[s]; + } +}; + +template +class MGPreconditioner : public LinearFunction< Lattice > { +public: + using LinearFunction >::operator(); + + typedef Aggregation Aggregates; + typedef typename Aggregation::FineField FineField; + typedef typename Aggregation::CoarseVector CoarseVector; + typedef typename Aggregation::CoarseMatrix CoarseMatrix; + typedef LinearOperatorBase FineOperator; + typedef LinearFunction FineSmoother; + typedef LinearOperatorBase CoarseOperator; + typedef LinearFunction CoarseSolver; + + Aggregates & _Aggregates; + FineOperator & _FineOperator; + FineSmoother & _PreSmoother; + FineSmoother & _PostSmoother; + CoarseOperator & _CoarseOperator; + CoarseSolver & _CoarseSolve; + CoarseSolver & _CoarseGuesser; + + int level; void Level(int lv) {level = lv; }; + + MGPreconditioner(Aggregates &Agg, + FineOperator &Fine, + FineSmoother &PreSmoother, + FineSmoother &PostSmoother, + CoarseOperator &CoarseOperator_, + CoarseSolver &CoarseSolve_, + CoarseSolver &CoarseGuesser_) + : _Aggregates(Agg), + _FineOperator(Fine), + _PreSmoother(PreSmoother), + _PostSmoother(PostSmoother), + _CoarseOperator(CoarseOperator_), + _CoarseSolve(CoarseSolve_), + _CoarseGuesser(CoarseGuesser_), + level(1) { } + + virtual void operator()(const FineField &in, FineField & out) + { + GridBase *CoarseGrid = _Aggregates.CoarseGrid; + CoarseVector Csrc(CoarseGrid); + CoarseVector Csol(CoarseGrid); + FineField vec1(in.Grid()); + FineField vec2(in.Grid()); + + double t; + out = Zero(); + t=-usecond(); + _PreSmoother(in,out); + t+=usecond(); + std::cout< +class ShiftedLinearOperator : public LinearOperatorBase { + LinearOperatorBase &_Op; + RealD shift; +public: + ShiftedLinearOperator(RealD _shift, LinearOperatorBase &Op) : shift(_shift), _Op(Op) {} + 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); } +}; + +template +void runMG( + GridCartesian *FGrid, + GridCartesian *Coarse5d, + GridCartesian *CoarseCoarse5d, + GridCartesian *CoarseCoarseCoarse5d, + GridCartesian *CoarseCoarseCoarseCoarse5d, + NextToNearestStencilGeometry5D geom, + PVdagM_t &PVdagM, + ShiftedPVdagM_t &ShiftedPVdagM, + Subspace &AggregatesPD +) { + std::vector subspace = AggregatesPD.subspace; + assert((int)subspace.size() == NB); + const int nbasis = NB; + const int cb = 0; + + CoarseVector c_src(Coarse5d); + CoarseVector c_res(Coarse5d); + Complex one(1.0); + + LatticeFermionD f_src(FGrid); + LatticeFermionD f_res(FGrid); + + TrivialPrecon simpleC; + TrivialPrecon simple_fine; + + ////////////////////////////////////////////////////////////////////// + // Level 0→1: coarsen PVdagM, build LinOpCoarse + ////////////////////////////////////////////////////////////////////// + LittleDiracOperator LittleDiracOpPV(geom, FGrid, Coarse5d); + LittleDiracOpPV.CoarsenOperator(PVdagM, AggregatesPD); + + NonHermitianLinearOperator LinOpCoarse(LittleDiracOpPV); + + ////////////////////////////////////////////////////////////////////// + // Baseline: plain PGCR on LinOpCoarse (reference for comparison) + ////////////////////////////////////////////////////////////////////// + std::cout< L2PGCR_baseline(3.0e-2,1100,LinOpCoarse,simpleC,10,10); + L2PGCR_baseline.Level(2); + L2PGCR_baseline.Name("Cbaseline"); + c_src = one; + c_res = Zero(); + L2PGCR_baseline(c_src,c_res); + + ////////////////////////////////////////////////////////////////////// + // psi_coarse: coarse projections of pre-GS fine null vectors. + // These are the Level 1 near-null vectors, promoted from Level 0. + // Used as the aggregation basis for Level 1→2 coarsening. + ////////////////////////////////////////////////////////////////////// + std::vector psi_coarse(nbasis, Coarse5d); + for (int k = 0; k < nbasis; k++) + AggregatesPD.ProjectToSubspace(psi_coarse[k], subspace[k]); + + ////////////////////////////////////////////////////////////////////// + // Diagnostics: W (fine projected matrix) and C (Galerkin check) + ////////////////////////////////////////////////////////////////////// + { + Eigen::MatrixXcd W = Eigen::MatrixXcd::Zero(nbasis, nbasis); + LatticeFermion ftmp(FGrid); + for (int j = 0; j < nbasis; j++) { + PVdagM.Op(subspace[j], ftmp); + for (int i = 0; i < nbasis; i++) + W(i,j) = TensorRemove(innerProduct(subspace[i], ftmp)); + } + RealD normW = W.norm(); + std::cout << GridLogMessage << "Fine projected matrix ||W|| = " << normW << std::endl; + + Eigen::MatrixXcd C = Eigen::MatrixXcd::Zero(nbasis, nbasis); + CoarseVector Ac(Coarse5d); + for (int l = 0; l < nbasis; l++) { + LinOpCoarse.Op(psi_coarse[l], Ac); + for (int k = 0; k < nbasis; k++) + C(k,l) = TensorRemove(innerProduct(psi_coarse[k], Ac)); + } + RealD normC = C.norm(); + RealD normCmCdag = (C - C.adjoint()).norm(); + std::cout << GridLogMessage << "Coarse null matrix ||C|| = " << normC << std::endl; + std::cout << GridLogMessage << "Coarse null matrix ||C - C†||/||C|| = " << normCmCdag/normC << std::endl; + std::cout << GridLogMessage << "Galerkin check ||C||/||W|| = " << normC/normW << std::endl; + } + + ////////////////////////////////////////////////////////////////////// + // Level 1→2: set up aggregation using psi_coarse as subspace. + // Block factor 2,2,3,2 (removes odd local sublattice in z given MPI + // geometry 3×6×4×4 where z-local at Level 1 is 6). + // psi_coarse are assigned directly; CoarsenOperator performs + // block-GS orthogonalisation before building LinOpCoarseCoarse. + ////////////////////////////////////////////////////////////////////// + // innerProduct(CoarseSiteObj, CoarseSiteObj) returns iScalar, so CComplex + // for the L1→L2 level must be iScalar, not vTComplex. + typedef typename CoarseVector::vector_object CoarseSiteObj; + typedef iScalar vTTComplex; + typedef GeneralCoarsenedMatrix LittleDiracOperatorL2; + typedef typename LittleDiracOperatorL2::CoarseVector CoarseCoarseVector; + typedef Aggregation SubspaceL2; + typedef MGPreconditioner L1to2MG; + + SubspaceL2 AggregatesL2(CoarseCoarse5d, Coarse5d, cb); + for (int k = 0; k < nbasis; k++) + AggregatesL2.subspace[k] = psi_coarse[k]; + + NextToNearestStencilGeometry5D geom2(CoarseCoarse5d); + LittleDiracOperatorL2 LittleDiracOpL2(geom2, Coarse5d, CoarseCoarse5d); + LittleDiracOpL2.CoarsenOperator(LinOpCoarse, AggregatesL2); + + NonHermitianLinearOperator LinOpCC(LittleDiracOpL2); + + TrivialPrecon simpleCC; + + ////////////////////////////////////////////////////////////////////// + // Lüscher deflation guesser for L3PGCR. + // Step 1: project psi_coarse[k] (promoted fine null vectors) to + // CoarseCoarseVector space — these cover the zero-momentum + // component of the near-null space of LinOpCC. + // Step 2: breed Nextra additional null vectors directly on LinOpCC + // using GCR with random sources — these pick up near-null + // modes at all spatial frequencies not spanned by step 1. + // Step 3: build C_{st} = over the + // full augmented basis and invert directly via Eigen LU. + ////////////////////////////////////////////////////////////////////// + std::vector psi_cc(nbasis, CoarseCoarse5d); + for (int k = 0; k < nbasis; k++) + AggregatesL2.ProjectToSubspace(psi_cc[k], psi_coarse[k]); + + { + int Nextra = nbasis; // breed as many extra as we have promoted ones + if ( getenv("CC_NEXTRA") ) Nextra = atoi(getenv("CC_NEXTRA")); + GridParallelRNG RNG_CC(CoarseCoarse5d); + RNG_CC.SeedFixedIntegers({11,13,17,19}); + PrecGeneralisedConjugateResidualNonHermitian + nullGCR(1e-2, 200, LinOpCC, simpleCC, 32, 32); + CoarseCoarseVector tmp(CoarseCoarse5d); + for (int k = 0; k < Nextra; k++) { + CoarseCoarseVector src(CoarseCoarse5d); + gaussian(RNG_CC, src); + tmp = Zero(); + nullGCR(src, tmp); + psi_cc.push_back(tmp); + } + std::cout << GridLogMessage << "LinOpCC deflation basis: " << nbasis + << " promoted + " << Nextra << " bred = " << psi_cc.size() << " total" << std::endl; + } + + const int Naug = psi_cc.size(); + Eigen::MatrixXcd Ccc = Eigen::MatrixXcd::Zero(Naug, Naug); + { + CoarseCoarseVector Acc(CoarseCoarse5d); + for (int l = 0; l < Naug; l++) { + LinOpCC.Op(psi_cc[l], Acc); + for (int k = 0; k < Naug; k++) + Ccc(k,l) = TensorRemove(innerProduct(psi_cc[k], Acc)); + } + } + { + RealD normCcc = Ccc.norm(); + RealD normCccmCdag = (Ccc - Ccc.adjoint()).norm(); + std::cout << GridLogMessage << "Coarse-coarse deflation matrix ||Ccc|| = " << normCcc << std::endl; + std::cout << GridLogMessage << "Coarse-coarse deflation matrix ||Ccc-Ccc†||/||Ccc|| = " << normCccmCdag/normCcc << std::endl; + } + Eigen::MatrixXcd Ccc_inv = Ccc.inverse(); + LuscherGuesser CCDeflGuesser(psi_cc, Ccc_inv); + + ////////////////////////////////////////////////////////////////////// + // Level 2→3: coarsen LinOpCC using the RAW promoted psi_cc as aggregation + // to build the Level 4 (coarse-coarse-coarse) operator. + // psi_cc[0..nbasis-1] are the coarse-coarse near-null vectors, projected + // from the RAW psi_coarse (themselves projected from the RAW fine null + // vectors) -- the pre-block-GS chain the whole construction depends on. + // CoarsenOperator block-GS orthogonalises AggregatesL3.subspace IN PLACE, + // so assign COPIES of psi_cc and keep psi_cc itself raw. + // + // Tensor depth deepens once more: innerProduct(CoarseCoarseSiteObj,...) returns + // iScalar, so CComplex for the L2→L3 level is iScalar>. + ////////////////////////////////////////////////////////////////////// + typedef typename CoarseCoarseVector::vector_object CoarseCoarseSiteObj; + typedef iScalar vTTTComplex; + typedef GeneralCoarsenedMatrix LittleDiracOperatorL3; + typedef typename LittleDiracOperatorL3::CoarseVector CoarseCoarseCoarseVector; + typedef Aggregation SubspaceL3; + typedef MGPreconditioner L2to3MG; + + SubspaceL3 AggregatesL3(CoarseCoarseCoarse5d, CoarseCoarse5d, cb); + for (int k = 0; k < nbasis; k++) + AggregatesL3.subspace[k] = psi_cc[k]; // raw promoted; COPY, keeps psi_cc raw + + NextToNearestStencilGeometry5D geom3(CoarseCoarseCoarse5d); + LittleDiracOperatorL3 LittleDiracOpL3(geom3, CoarseCoarse5d, CoarseCoarseCoarse5d); + LittleDiracOpL3.CoarsenOperator(LinOpCC, AggregatesL3); // block-GS's AggregatesL3.subspace in place + + NonHermitianLinearOperator LinOpCCC(LittleDiracOpL3); + TrivialPrecon simpleCCC; + + ////////////////////////////////////////////////////////////////////// + // Level 3→4: coarsen LinOpCCC to build the Level 5 operator, using a + // TRUNCATED basis of only the first NB5 (< nbasis) raw promoted null vectors. + // psi_ccc[k] = raw psi_cc projected through the (block-GS'd) L3 aggregation + // -- the pre-block-GS chain continued one level deeper. We keep only the + // leading NB5: after the global orthogonalisation of the original fine null + // vectors the early indices retain the most-null content (shared low-mode + // components are peeled in first), so the leading NB5 are the crudely-most- + // null slice. This is the cheap "first 30" truncation test; a principled + // sigma-ordered rotation of psi_ccc would replace the slice, not the idea. + // NB: a positive result is conservative (sigma-ordering can only help); a + // negative one is inconclusive until the sigma-ordered NB5 is tried. + // + // Tensor depth deepens once more: CComplex for the L3→L4 level is + // iScalar. NB5 (the coarse dimension) is independent of the + // depth -- it just makes the coarsest site vector NB5-dimensional. + ////////////////////////////////////////////////////////////////////// + const int NB5 = 30; // compile-time: changing it re-instantiates the L4/L5 tensors + std::cout << GridLogMessage << "PARAM NB5 (truncated coarsest basis) = " << NB5 << std::endl; + assert(NB5 <= nbasis); + + std::vector psi_ccc(nbasis, CoarseCoarseCoarse5d); + for (int k = 0; k < nbasis; k++) + AggregatesL3.ProjectToSubspace(psi_ccc[k], psi_cc[k]); // raw psi_cc -> L4 null vectors + + ////////////////////////////////////////////////////////////////////// + // Optional sigma-ordering of psi_ccc (SVD_REORDER set): replace the crude + // first-NB5 slice with the NB5 genuinely-most-null directions of span(psi_ccc) + // under LinOpCCC. For a NON-NORMAL operator the nullness measure is the + // singular value of A restricted to the span -- eig of Q†A†AQ -- NOT the + // numerical range Q†AQ (which non-normality contaminates). Robust route: + // whiten by the Gram (drop near-dependent directions), Hermitian-eig the + // whitened A†A, rotate. The printed singular spectrum IS the SVD study: where + // it falls off tells you the natural NB5, and the same numbers illuminate why + // the earlier singular-subspace deflation re-entered. Safe here because we + // ORDER vectors that then feed a Galerkin projection, not REMOVE a subspace. + // Default (unset) leaves psi_ccc in raw order == the "first 30" test. + ////////////////////////////////////////////////////////////////////// + if ( getenv("SVD_REORDER") ) { + std::cout << GridLogMessage << "SVD_REORDER: sigma-ordering psi_ccc under LinOpCCC" << std::endl; + + Eigen::MatrixXcd G(nbasis,nbasis); // Gram = Psi^dag Psi + for (int i=0;i Apsi(nbasis, CoarseCoarseCoarse5d); + for (int j=0;j tol*max; T = Ug diag(1/sqrt g). + // Q = Psi T is then orthonormal (Q^dag Q = T^dag G T = I). + Eigen::SelfAdjointEigenSolver esG(G); + Eigen::VectorXd g = esG.eigenvalues(); // ascending, real + RealD gmax = g(nbasis-1); + RealD gtol = 1.0e-9 * gmax; + int keep = 0; for (int i=0;i gtol) keep++; + std::cout << GridLogMessage << " Gram spectrum: min=" << g(0) << " max=" << gmax + << " cond=" << gmax/std::max(g(0),1.0e-300) << " keep=" << keep << "/" << nbasis << std::endl; + assert(keep >= NB5); + + Eigen::MatrixXcd T(nbasis, keep); // whitening (largest-g first) + { int c=0; + for (int i=nbasis-1;i>=0;i--) if (g(i) > gtol) { T.col(c) = esG.eigenvectors().col(i)/std::sqrt(g(i)); c++; } + } + + Eigen::MatrixXcd Mw = T.adjoint() * M * T; // whitened A^dagA (keep x keep, Hermitian) + Eigen::SelfAdjointEigenSolver esM(Mw); + Eigen::VectorXd s2 = esM.eigenvalues(); // ascending sigma^2 (most-null first) + std::cout << GridLogMessage << " Singular spectrum sigma_k (most-null first):" << std::endl; + for (int k=0;k vTTTTComplex; + typedef GeneralCoarsenedMatrix LittleDiracOperatorL4; + typedef typename LittleDiracOperatorL4::CoarseVector CoarseCoarseCoarseCoarseVector; + typedef Aggregation SubspaceL4; + typedef MGPreconditioner L3to4MG; + + SubspaceL4 AggregatesL4(CoarseCoarseCoarseCoarse5d, CoarseCoarseCoarse5d, cb); + for (int k = 0; k < NB5; k++) + AggregatesL4.subspace[k] = psi_ccc[k]; // FIRST NB5 raw promoted vectors (truncation) + + NextToNearestStencilGeometry5D geom4(CoarseCoarseCoarseCoarse5d); + LittleDiracOperatorL4 LittleDiracOpL4(geom4, CoarseCoarseCoarse5d, CoarseCoarseCoarseCoarse5d); + LittleDiracOpL4.CoarsenOperator(LinOpCCC, AggregatesL4); // block-GS's AggregatesL4.subspace in place + + NonHermitianLinearOperator LinOpCCCC(LittleDiracOpL4); + TrivialPrecon simpleCCCC; + + ////////////////////////////////////////////////////////////////////// + // Level 5 bottom solve: GCR on a SHIFTED LinOpCCCC (the coarsest, most + // non-normal operator). l5_shift slides its field of values off the origin; + // defaults to 0.0 (bare LinOpCCCC) until opted in. This is the level a dense + // direct inverse would eventually replace: rank = NB5 * sites(clatt4). + ////////////////////////////////////////////////////////////////////// + RealD l5_shift = 0.0; + if(getenv("l5_shift")) l5_shift = atof(getenv("l5_shift")); + std::cout << GridLogMessage << "PARAM l5_shift = " << l5_shift << std::endl; + + ShiftedLinearOperator ShiftedLinOpCCCC(l5_shift, LinOpCCCC); + PrecGeneralisedConjugateResidualNonHermitian L5PGCR(1.0e-1,200,ShiftedLinOpCCCC,simpleCCCC,16,16); + L5PGCR.Level(5); + L5PGCR.Name("CCCCouter"); + + ////////////////////////////////////////////////////////////////////// + // Level 3→4 V-cycle: depth-2 SHIFTED smoother on LinOpCCC + Level 5 bottom. + // Level 4 is no longer the bottom -- it is smoothed shallowly and recursed to + // Level 5, mirroring how Level 3 recurses to Level 4. + ////////////////////////////////////////////////////////////////////// + RealD ccc_smoother_shift = 0.05; + int ccc_smoother_nstep = 2; + if(getenv("ccc_smoother_shift")) ccc_smoother_shift = atof(getenv("ccc_smoother_shift")); + if(getenv("ccc_smoother_nstep")) ccc_smoother_nstep = atoi(getenv("ccc_smoother_nstep")); + + ShiftedLinearOperator ShiftedLinOpCCC(ccc_smoother_shift, LinOpCCC); + PrecGeneralisedConjugateResidualNonHermitian + CoarseCoarseCoarseSmootherGCR(0.01,1,ShiftedLinOpCCC,simpleCCC,ccc_smoother_nstep,ccc_smoother_nstep); + CoarseCoarseCoarseSmootherGCR.SetZeroGuess(1); // smoother slot: caller zeroes guess + CoarseCoarseCoarseSmootherGCR.Level(4); + CoarseCoarseCoarseSmootherGCR.Name("CCCsmoother"); + + L3to4MG L3to4Precon(AggregatesL4, + LinOpCCC, + simpleCCC, // no pre-smoother + CoarseCoarseCoarseSmootherGCR, // post-smoother: depth-2 shifted GCR + LinOpCCCC, + L5PGCR, + simpleCCCC); // trivial guesser at the bottom + + ////////////////////////////////////////////////////////////////////// + // Level 4 (coarse-coarse-coarse) solve: GCR preconditioned by the L3→L4 V-cycle. + ////////////////////////////////////////////////////////////////////// + PrecGeneralisedConjugateResidualNonHermitian L4MGsolver(1.0e-1,200,LinOpCCC,L3to4Precon,16,16); + L4MGsolver.Level(4); + L4MGsolver.Name("CCCouter"); + + ////////////////////////////////////////////////////////////////////// + // Level 2→3 V-cycle: depth-2 SHIFTED smoother on LinOpCC + Level 4 solve. + // The shift slides the coarse-coarse field of values off the origin so a + // 2-step smoother has something to bite on a non-normal operator (IRS idea). + ////////////////////////////////////////////////////////////////////// + RealD cc_smoother_shift = 0.01; + int cc_smoother_nstep = 2; + if(getenv("cc_smoother_shift")) cc_smoother_shift = atof(getenv("cc_smoother_shift")); + if(getenv("cc_smoother_nstep")) cc_smoother_nstep = atoi(getenv("cc_smoother_nstep")); + + ShiftedLinearOperator ShiftedLinOpCC(cc_smoother_shift, LinOpCC); + PrecGeneralisedConjugateResidualNonHermitian + CoarseCoarseSmootherGCR(0.01,1,ShiftedLinOpCC,simpleCC,cc_smoother_nstep,cc_smoother_nstep); + CoarseCoarseSmootherGCR.SetZeroGuess(1); // smoother slot: caller zeroes guess + CoarseCoarseSmootherGCR.Level(3); + CoarseCoarseSmootherGCR.Name("CCsmoother"); + + L2to3MG L2to3Precon(AggregatesL3, + LinOpCC, + simpleCC, // no pre-smoother + CoarseCoarseSmootherGCR, // post-smoother: depth-2 shifted GCR + LinOpCCC, + L4MGsolver, // coarse solve is now the L3→L4 V-cycle + simpleCCC); // trivial guesser + + ////////////////////////////////////////////////////////////////////// + // Level 3 (coarse-coarse) solve: GCR preconditioned by the L2→L3 V-cycle. + // Replaces the plain L3PGCR of the 3-level build -- the coarse-coarse level + // is now smoothed shallowly and recursed rather than solved deeply. + ////////////////////////////////////////////////////////////////////// + PrecGeneralisedConjugateResidualNonHermitian L3MGsolver(1.0e-1,200,LinOpCC,L2to3Precon,16,16); + L3MGsolver.Level(3); + L3MGsolver.Name("CCouter"); + + ////////////////////////////////////////////////////////////////////// + // Coarse-level GCR smoother for Level 1→2 V-cycle. + // Mirrors fine-grid SmootherGCR: shifted operator + fixed step count. + // coarse_smoother_shift and coarse_smoother_nstep are the tuning knobs. + ////////////////////////////////////////////////////////////////////// + RealD coarse_smoother_shift = 0.01; + int coarse_smoother_nstep = 2; // depth-2 smoother on the coarse level + if(getenv("coarse_smoother_shift")) coarse_smoother_shift = atof(getenv("coarse_smoother_shift")); + if(getenv("coarse_smoother_nstep")) coarse_smoother_nstep = atoi(getenv("coarse_smoother_nstep")); + + ShiftedLinearOperator ShiftedLinOpCoarse(coarse_smoother_shift, LinOpCoarse); + PrecGeneralisedConjugateResidualNonHermitian CoarseSmootherGCR(0.01,1,ShiftedLinOpCoarse,simpleC,coarse_smoother_nstep,coarse_smoother_nstep); + CoarseSmootherGCR.SetZeroGuess(1); // smoother slot: caller zeroes guess + CoarseSmootherGCR.Level(2); + CoarseSmootherGCR.Name("Csmoother"); + + ////////////////////////////////////////////////////////////////////// + // Level 1→2 V-cycle preconditioner. + ////////////////////////////////////////////////////////////////////// + L1to2MG L1to2Precon(AggregatesL2, + LinOpCoarse, + simpleC, // no pre-smoother (matches fine-grid setup) + CoarseSmootherGCR, // post-smoother: depth-2 shifted GCR + LinOpCC, + L3MGsolver, // coarse-coarse solve is now the L2→L3 V-cycle + CCDeflGuesser); // Lüscher guesser: psi_cc C^{-1} psi_cc† + + ////////////////////////////////////////////////////////////////////// + // Standalone Level 1 two-level solve test. + // Compare against plain PGCR baseline above. + ////////////////////////////////////////////////////////////////////// + std::cout< L2MGsolver(3.0e-2,200,LinOpCoarse,L1to2Precon,16,16); + L2MGsolver.Level(2); + L2MGsolver.Name("Couter"); + c_res = Zero(); + L2MGsolver(c_src,c_res); + + std::cout << GridLogMessage << "Level 1 two-level test: PVdagM operator uses:" << std::endl; + PVdagM.getApplications(); + PVdagM.clear(); + + ////////////////////////////////////////////////////////////////////// + // Full five-level outer solve + ////////////////////////////////////////////////////////////////////// + std::cout< SmootherGCR(0.01,1,ShiftedPVdagM,simple_fine,16,16); + SmootherGCR.SetZeroGuess(1); // pre/post smoother slots zero their guess + SmootherGCR.Level(1); + SmootherGCR.Name("Fsmoother"); + + f_src = one; + + // Pre-smoother: none (TrivialPrecon); post-smoother: shifted PGCR. + // Coarse solver: L2MGsolver (PGCR preconditioned by Level 1→2 V-cycle). + TwoLevelMG ThreeLevelPrecon(AggregatesPD, + PVdagM, + simple_fine, + SmootherGCR, + LinOpCoarse, + L2MGsolver, + simpleC); + + PrecGeneralisedConjugateResidualNonHermitian L1PGCR(1.0e-8,1000,PVdagM,ThreeLevelPrecon,16,16); + L1PGCR.Level(1); + L1PGCR.Name("Fouter"); + + f_res = Zero(); + L1PGCR(f_src,f_res); + + std::cout << GridLogMessage << "Five-level outer solve: PVdagM operator uses:" << std::endl; + PVdagM.getApplications(); + PVdagM.clear(); +} + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + const int Ls = 24; + RealD M5 = 1.8; + RealD b = 1.5; + RealD c = 0.5; + RealD mass = 0.00078; + if ( getenv("MASS") ) mass = atof(getenv("MASS")); + + const int nbasis = 60; + + std::cout << GridLogMessage << "Mass: " << mass << ", Ls: " << Ls << ", b=" << b << ", c=" << c << std::endl; + std::cout << GridLogMessage << "nbasis: " << nbasis << std::endl; + + 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); + + // Level 1 coarse grid: block 2^4 from fine (48×48×48×96 → 24×24×24×48, Ls=1) + Coordinate clatt = lat_size; + for (int d = 0; d < 4; d++) clatt[d] /= 2; + std::cout << GridLogMessage << "Level 1 coarse lattice: " << clatt << std::endl; + + GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); + GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d); + + // Level 2 coarse-coarse grid: block 2,2,3,3 from Level 1 (24×24×24×48 → 12×12×8×16, Ls=1). + // MPI geometry 3.6.4.4 (288 ranks): fine local {16,8,12,24}. + // Level 1 local {8,4,6,12}; Level 2 local {4,2,2,4}. + // z blocked by 3: z-Level1-local=6; 6/3=2 (even), 6/2=3 (odd) → must use 3. + // t blocked by 3: t-Level1-local=12; 12/3=4 divisible by Nsimd=4 (gen-simd-width=64). + // t-block=2 gives t2-local=6, 6 mod 4 ≠ 0, fails Grid SIMD assertion. ✓ + // With {4,2,2,4}: Nsimd=4 goes into x or t (both =4). ✓ + Coordinate clatt2 = clatt; + clatt2[0] /= 2; + clatt2[1] /= 2; + clatt2[2] /= 3; + clatt2[3] /= 3; + std::cout << GridLogMessage << "Level 2 coarse-coarse lattice: " << clatt2 << std::endl; + + GridCartesian *CoarseCoarse4d = SpaceTimeGrid::makeFourDimGrid(clatt2, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); + GridCartesian *CoarseCoarse5d = SpaceTimeGrid::makeFiveDimGrid(1,CoarseCoarse4d); + + // Level 3 coarse-coarse-coarse grid: block clatt2 = {12,12,8,16} -> {6,12,8,8}. + // GEOMETRY (mpi 3.6.4.4, Nsimd=4 => SIMD layout {1,1,2,2}, factor 2 on z and t): + // every grid needs z-local and t-local EVEN. clatt2-local is {4,2,2,4}, so + // z-local=2 is already at its minimum even value and CANNOT be blocked (2->1 + // is odd and trips the SIMD assertion); y-local=2 would go to 1 (degenerate). + // Only x and t have room, so block {2,1,1,2}: clatt3 {6,12,8,8}, L4-local + // {2,2,2,2} -- all dims even and >=2. z stays unblocked by construction. + Coordinate clatt3 = clatt2; + clatt3[0] /= 2; // x: 12 -> 6 (x-local 4 -> 2) + // clatt3[1] (y) unblocked: y-local 2 -> blocking gives 1 (degenerate) + // clatt3[2] (z) unblocked: z-local 2 is SIMD-pinned even, cannot halve + clatt3[3] /= 2; // t: 16 -> 8 (t-local 4 -> 2) + std::cout << GridLogMessage << "Level 3 coarse-coarse-coarse lattice: " << clatt3 << std::endl; + + GridCartesian *CoarseCoarseCoarse4d = SpaceTimeGrid::makeFourDimGrid(clatt3, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); + GridCartesian *CoarseCoarseCoarse5d = SpaceTimeGrid::makeFiveDimGrid(1,CoarseCoarseCoarse4d); + + // Level 4 coarse^4 grid: block clatt3 = {6,12,8,8} -> {3,6,8,8} via {2,2,1,1}. + // mpi 3.6.4.4 => clatt4-local {1,1,2,2}: z-local=2, t-local=2 stay EVEN (SIMD + // factor 2 pins them), so z,t are unblocked; x,y (SIMD factor 1) halve to + // local 1 -- fully distributed but legal for the halo-depth-1 NextToNearest + // stencil. 1152 sites; with NB5=30 that is the 34,560-rank coarsest operator + // a dense direct inverse would target. + Coordinate clatt4 = clatt3; + clatt4[0] /= 2; // x: 6 -> 3 (x-local 2 -> 1) + clatt4[1] /= 2; // y: 12 -> 6 (y-local 2 -> 1) + // clatt4[2] (z) unblocked: z-local 2 is SIMD-pinned even + // clatt4[3] (t) unblocked: t-local 2 is SIMD-pinned even + std::cout << GridLogMessage << "Level 4 coarse^4 lattice: " << clatt4 << std::endl; + + GridCartesian *CoarseCoarseCoarseCoarse4d = SpaceTimeGrid::makeFourDimGrid(clatt4, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); + GridCartesian *CoarseCoarseCoarseCoarse5d = SpaceTimeGrid::makeFiveDimGrid(1,CoarseCoarseCoarseCoarse4d); + + std::vector seeds4({1,2,3,4}); + std::vector seeds5({5,6,7,8}); + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + + LatticeGaugeField Umu(UGrid); + std::cout << GridLogMessage << "Reading gauge field" << std::endl; + FieldMetaData header; + std::string file("/ccs/home/poare/ckpoint_lat.1000"); + NerscIO::readConfiguration(Umu,header,file); + + RealD b_ = 1.5; + RealD c_ = 0.5; + 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; + typedef ShiftedPVdagMLinearOperator ShiftedPVdagM_t; + typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef LittleDiracOperator::CoarseVector CoarseVector; + typedef Aggregation Subspace; + typedef MGPreconditioner TwoLevelMG; + + PVdagM_t PVdagM(Ddwf,Dpv); + ShiftedPVdagM_t ShiftedPVdagM(0.01,Ddwf,Dpv); + + NextToNearestStencilGeometry5D geom(Coarse5d); + + // Subspace cache: save after generation, reload on subsequent runs to skip expensive setup. + // Set SUBSPACE_FILE to override the default path. + 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")); + + // Check if subspace file exists (boss rank checks, result broadcast via GlobalSum). + 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); + // Insurance: GLOBAL (whole-lattice) orthonormalise, in case the cached file + // predates the GlobalOrthonormalise() that CreateSubspaceGCR now applies + // (Aggregates.h:196). It is span-preserving and makes the vectors globally + // orthonormal -- it is NOT the block Orthogonalise() below, so it does NOT + // cause the psi_coarse->e_k trap. It also (re)establishes the weak nullness + // gradient (shared most-null components peeled into the early indices) that + // the "first NB5" truncation relies on. Idempotent if the file was already + // globally orthonormal. The RAW subspace copy in runMG happens AFTER this + // call, so the raw-null (pre-block-GS) discipline is preserved. + AggregatesGCR.GlobalOrthonormalise(); + // DO NOT block-orthogonalise here: runMG copies subspace[] as the RAW + // (pre-block-GS) basis and CoarsenOperator block-GS's it in place later. + // Orthogonalising now defeats the raw-null discipline (psi_coarse -> e_k) + // and poisons L2/L3/L4. See project_block_orthogonalise_leak. + // AggregatesGCR.Orthogonalise(); + std::cout << GridLogMessage << "Subspace loaded, globally orthonormalised (raw block basis preserved)." << std::endl; + } else { + std::cout << GridLogMessage << "*** GCR subspace generation ***" << std::endl; + AggregatesGCR.CreateSubspaceGCR(RNG5,PVdagM,nbasis); + std::cout << GridLogMessage << "Subspace generation: PVdagM operator uses:" << std::endl; + PVdagM.getApplications(); + PVdagM.clear(); + saveSubspace(AggregatesGCR.subspace, subspace_file); + std::cout << GridLogMessage << "Subspace saved to: " << subspace_file << std::endl; + } + + runMG( + FGrid, + Coarse5d, + CoarseCoarse5d, + CoarseCoarseCoarse5d, + CoarseCoarseCoarseCoarse5d, + geom, + PVdagM, + ShiftedPVdagM, + AggregatesGCR + ); + + std::cout << GridLogMessage << "Done" << std::endl; + Grid_finalize(); + return 0; +} diff --git a/examples/Example_pvdagm_census.cc b/examples/Example_pvdagm_census.cc new file mode 100644 index 000000000..e04b81ed7 --- /dev/null +++ b/examples/Example_pvdagm_census.cc @@ -0,0 +1,717 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./examples/Example_pvdagm_census.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 */ + +// Spectral census of the coarsened PVdagM operator A_c. +// +// Measures the three sets that discriminate between the candidate explanations +// for slow coarse-grid Krylov convergence: +// +// 0. Fine Ritz diagonal of RAW subspace vectors (pre-block-orthog). +// NB CoarsenOperator block-orthogonalises subspace[] IN PLACE; all +// nulliness/deflation bases must be built from a raw copy. +// 1. Adjoint correctness check == (fail-fast) +// 2. Raw-vector coarse images vs A_c: RQ (must equal CENSUS 0 by Galerkin), +// ||A_c psi_c||/||psi_c||, and representability error +// 3. sigma_max^2 = lambda_max(A_c^dag A_c) via power method +// 4. Low singular values Chebyshev-filtered IRL on A_c^dag A_c +// -> sigma_min census = pseudospectrum of A_c evaluated at the origin +// 5. Half-plane margin lambda_min/max of H = (A_c + A_c^dag)/2 +// -> min Re W(A_c); positive-real check (Eisenstat-Elman-Schultz bound) +// +// Interpretation: +// sigma_min ~ min|lambda|, ~nbasis tiny then gap : effectively normal, bipartite +// sigma_min ~ min|lambda|, dense low tail : normal but rank-starved +// sigma_min << min|lambda| : non-normal near origin +// lambda_min(H) < 0 : half-plane condition violated +// +// Requires the dagger code path in GeneralCoarsenedMatrix: +// _Adag allocated, PopulateAdag active, _Adag exchanged, hermitian=0. +// +// Env vars: +// MASS fermion mass (default 0.00078) +// SUBSPACE_FILE subspace cache path +// CoarseSolverShift shift baked into coarsening (default 0.0: pure Galerkin) +// CENSUS_NSTOP converged low modes wanted (default 60) +// CENSUS_NK Lanczos Nk (default 96) +// CENSUS_NM Lanczos Nm (default 192) +// CENSUS_TOL Lanczos residual (default 1e-5) +// CENSUS_MAXIT Lanczos max restarts (default 50) +// CHEBY_LO filter low edge in sigma^2 (default 4.0) +// CHEBY_HI filter high edge; 0 = auto from power method x1.1 +// CHEBY_ORDER filter order (default 401) +// filter gain at 0 ~ cosh(order*2*sqrt(lo/hi)); with +// hi~2200, lo=4, order=401 => gain ~ 1e14. lo=0.01 at +// order 201 gives gain ~1.4 (stagnation). + +#include +#include +#include +#include + +using namespace std; +using namespace Grid; + +RealD mass = 0.00078; +RealD CoarseSolverShift = 0.0; +int CensusNstop = 60; +int CensusNk = 96; +int CensusNm = 192; +RealD CensusTol = 1.0e-5; +int CensusMaxIt = 50; +RealD ChebyLo = 4.0; // sigma^2 cutoff: amplifies sigma < 2. Filter gain ~ cosh(order*2*sqrt(lo/hi)) +RealD ChebyHi = 0.0; // 0 => auto: 1.1 * power-method sigma_max^2 +int ChebyOrder = 401; +RealD CGdeflTol = 1.0e-8; // CENSUS 6 deflated-CG tolerance +int CGdeflMaxIt = 4000; // CENSUS 6 deflated-CG max iterations +int DeflRank = 0; // CENSUS 6 deflation rank; 0 => all available per basis + +void ParseEnvironment(void) +{ + if(getenv("MASS")) mass = atof(getenv("MASS")); + if(getenv("CoarseSolverShift")) CoarseSolverShift = atof(getenv("CoarseSolverShift")); + if(getenv("CENSUS_NSTOP")) CensusNstop = atoi(getenv("CENSUS_NSTOP")); + if(getenv("CENSUS_NK")) CensusNk = atoi(getenv("CENSUS_NK")); + if(getenv("CENSUS_NM")) CensusNm = atoi(getenv("CENSUS_NM")); + if(getenv("CENSUS_TOL")) CensusTol = atof(getenv("CENSUS_TOL")); + if(getenv("CENSUS_MAXIT")) CensusMaxIt = atoi(getenv("CENSUS_MAXIT")); + if(getenv("CHEBY_LO")) ChebyLo = atof(getenv("CHEBY_LO")); + if(getenv("CHEBY_HI")) ChebyHi = atof(getenv("CHEBY_HI")); + if(getenv("CHEBY_ORDER")) ChebyOrder = atoi(getenv("CHEBY_ORDER")); + if(getenv("CGDEFL_TOL")) CGdeflTol = atof(getenv("CGDEFL_TOL")); + if(getenv("CGDEFL_MAXIT")) CGdeflMaxIt = atoi(getenv("CGDEFL_MAXIT")); + if(getenv("DEFL_RANK")) DeflRank = atoi(getenv("DEFL_RANK")); + + std::cout << GridLogMessage << "PARAM: MASS " << mass << std::endl; + std::cout << GridLogMessage << "PARAM: CoarseSolverShift " << CoarseSolverShift << std::endl; + std::cout << GridLogMessage << "PARAM: CENSUS_NSTOP " << CensusNstop << std::endl; + std::cout << GridLogMessage << "PARAM: CENSUS_NK " << CensusNk << std::endl; + std::cout << GridLogMessage << "PARAM: CENSUS_NM " << CensusNm << std::endl; + std::cout << GridLogMessage << "PARAM: CENSUS_TOL " << CensusTol << std::endl; + std::cout << GridLogMessage << "PARAM: CENSUS_MAXIT " << CensusMaxIt << std::endl; + std::cout << GridLogMessage << "PARAM: CHEBY_LO " << ChebyLo << std::endl; + std::cout << GridLogMessage << "PARAM: CHEBY_HI " << ChebyHi << std::endl; + std::cout << GridLogMessage << "PARAM: CHEBY_ORDER " << ChebyOrder << 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); + } +}; + +// H = (A + A^dag)/2 : Hermitian part of the coarse operator. +// lambda_min(H) = min Re W(A) is the half-plane margin; the EES GCR +// convergence theorem requires it positive. +template +class HermitianPartOperator : public LinearOperatorBase { + Matrix &_Mat; +public: + HermitianPartOperator(Matrix &Mat): _Mat(Mat) {}; + 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){ HermOp(in,out); } + void AdjOp (const Field &in, Field &out){ HermOp(in,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()); + _Mat.M(in,out); + _Mat.Mdag(in,tmp); + out = 0.5*(out + tmp); + } +}; + +// s*I - Op : power method on this gives s - lambda_min(Op) for Hermitian Op. +template +class ShiftedNegatedOperator : public LinearOperatorBase { + LinearOperatorBase &_Op; + RealD s; +public: + ShiftedNegatedOperator(RealD _s, LinearOperatorBase &Op): _Op(Op), s(_s) {}; + 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){ HermOp(in,out); } + void AdjOp (const Field &in, Field &out){ HermOp(in,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){ + _Op.HermOp(in,out); + out = s*in - out; + } +}; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + ParseEnvironment(); + + const int Ls=24; + RealD M5=1.8; + RealD b=1.5; + RealD c=0.5; + const int nbasis = 60; + + std::cout << GridLogMessage << "Census of coarse PVdagM: mass=" << mass << " Ls=" << Ls << " nbasis=" << nbasis << std::endl; + + 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); + + // Blocking: default matches Example_pvdagm.cc; override with e.g. BLOCK=2.2.2.2 + Coordinate clatt = lat_size; + Coordinate Block({4,4,6,4}); + if ( getenv("BLOCK") ) { + GridCmdOptionIntVector(std::string(getenv("BLOCK")),Block); + GRID_ASSERT(Block.size()==4); + } + for(int d=0;d PVdagM_t; + typedef ShiftedPVdagMLinearOperator ShiftedPVdagM_t; + typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef LittleDiracOperator::CoarseVector CoarseVector; + typedef Aggregation Subspace; + + PVdagM_t PVdagM(Ddwf,Dpv); + ShiftedPVdagM_t ShiftedPVdagM(CoarseSolverShift,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); + } + + ////////////////////////////////////////////////////////////////////// + // Keep the RAW (pre-block-orthogonalisation) near-null vectors. + // CoarsenOperator block-orthogonalises subspace[] IN PLACE, after which + // subspace[k] is the orthonormal basis phi_k and Project(phi_k) = e_k, + // the block-constant unit vector -- NOT a near-null direction. + // All nulliness measurements and any deflation basis must use raw[]. + ////////////////////////////////////////////////////////////////////// + std::vector raw(nbasis,FGrid); + for(int k=0;k/ ~ the nulliness achieved at generation + // (~2e-3). O(0.1-10) values mean the cache holds orthogonalised vectors + // and must be regenerated. + ////////////////////////////////////////////////////////////////////// + std::cout << GridLogMessage << "=================================================" << std::endl; + std::cout << GridLogMessage << "CENSUS 0: fine Ritz diagonal of raw subspace vectors" << std::endl; + std::cout << GridLogMessage << "=================================================" << std::endl; + { + LatticeFermionD Ap(FGrid); + for(int k=0;k/ = " << rq + << " ||A psi||/||psi|| = " << std::sqrt(norm2(Ap)/n2psi) << std::endl; + } + } + + ////////////////////////////////////////////////////////////////////// + // Coarsen. hermitian=0 is REQUIRED: enables PopulateAdag so that + // Mdag applies A^dag rather than silently aliasing to A. + ////////////////////////////////////////////////////////////////////// + LittleDiracOperator LittleDiracOpPV(geom,FGrid,Coarse5d,0); + if ( CoarseSolverShift != 0.0 ) { + std::cout << GridLogMessage << "Coarsening SHIFTED operator, shift=" << CoarseSolverShift << std::endl; + LittleDiracOpPV.CoarsenOperator(ShiftedPVdagM, AggregatesGCR); + } else { + std::cout << GridLogMessage << "Coarsening pure Galerkin operator (no shift)" << std::endl; + LittleDiracOpPV.CoarsenOperator(PVdagM, AggregatesGCR); + } + + CoarseVector c_x(Coarse5d); + CoarseVector c_y(Coarse5d); + CoarseVector c_t1(Coarse5d); + CoarseVector c_t2(Coarse5d); + + ////////////////////////////////////////////////////////////////////// + // CENSUS 1: adjoint correctness (fail fast) + // == for random x,y + ////////////////////////////////////////////////////////////////////// + std::cout << GridLogMessage << "=================================================" << std::endl; + std::cout << GridLogMessage << "CENSUS 1: adjoint correctness of dagger code path" << std::endl; + std::cout << GridLogMessage << "=================================================" << std::endl; + random(CRNG,c_x); + random(CRNG,c_y); + LittleDiracOpPV.M (c_x,c_t1); // A x + LittleDiracOpPV.Mdag(c_y,c_t2); // A^dag y + ComplexD ip1 = innerProduct(c_y,c_t1); // + ComplexD ip2 = innerProduct(c_t2,c_x); // + RealD reldiff = abs(ip1-ip2)/abs(ip1); + std::cout << GridLogMessage << "CENSUS: = " << ip1 << std::endl; + std::cout << GridLogMessage << "CENSUS: = " << ip2 << std::endl; + std::cout << GridLogMessage << "CENSUS: rel diff = " << reldiff << " (expect ~1e-14; FAIL if O(1))" << std::endl; + GRID_ASSERT(reldiff < 1.0e-8); + + // Coarse near-null ("global") vectors psi_c[k] = P^dag raw[k], stored for the + // Ritz-matrix + deflation study in CENSUS 6 (filled in CENSUS 2's projection + // loop below, before raw[]/subspace[] are freed). + std::vector psi_c(nbasis,Coarse5d); + + ////////////////////////////////////////////////////////////////////// + // CENSUS 2: nulliness of the RAW vectors' coarse images against A_c. + // psi_c[k] = P^dag raw[k]. Galerkin guarantees the Rayleigh quotient + // equals CENSUS 0's fine value exactly (raw[k] is in span of its own + // chopped pieces) -- agreement is a machine-precision validation of + // the coarsening. ||A_c psi_c||/||psi_c|| is the sigma-relevant norm. + // The representability column ||raw - P psi_c||/||raw|| must be ~eps. + ////////////////////////////////////////////////////////////////////// + std::cout << GridLogMessage << "=================================================" << std::endl; + std::cout << GridLogMessage << "CENSUS 2: raw-vector coarse images against coarse operator" << std::endl; + std::cout << GridLogMessage << "=================================================" << std::endl; + { + LatticeFermionD back(FGrid); + for(int k=0;k/ = " << rq + << " ||A psi||/||psi|| = " << std::sqrt(n2Apsi/n2psi) + << " represent_err = " << represent << std::endl; + } + } + + // Fine subspace + raw copy are needed only through CENSUS 2; CENSUS 3-5 are + // entirely coarse (LittleDiracOpPV only), and the CENSUS 4 evec save writes the + // coarse vectors directly. Release the ~2*nbasis fine 5D fields (~14 GB/GCD at + // 2^4) HERE, before the order-ChebyOrder Lanczos whose padded coarse temporaries + // otherwise push host memory over the top on top of _A + _Adag (the AccCache + // CpuPtr!=NULL abort seen mid-iteration). + // Direct orthonormality check of the fine near-null vectors (GlobalOrthonormalise + // in CreateSubspaceGCR). raw is freed just below, so this runs here, not CENSUS 6. + // If this is ~0 but the coarse Gram S (CENSUS 6) is not, the gap is representability, + // not orthonormality. + { + Eigen::MatrixXcd Gfine(nbasis,nbasis); + for(int i=0;i(g.real(),g.imag()); + Gfine(j,i) = std::conj(Gfine(i,j)); + } + } + double GmI = (Gfine - Eigen::MatrixXcd::Identity(nbasis,nbasis)).norm(); + std::cout << GridLogMessage << "CENSUS 2b: fine Gram || - I||_F = " << GmI + << " (expect ~0 if fine vectors orthonormal)" << std::endl; + } + + raw.clear(); raw.shrink_to_fit(); + AggregatesGCR.subspace.clear(); AggregatesGCR.subspace.shrink_to_fit(); + + ////////////////////////////////////////////////////////////////////// + // CENSUS 3: sigma_max^2 = lambda_max( A_c^dag A_c ) by power method + ////////////////////////////////////////////////////////////////////// + std::cout << GridLogMessage << "=================================================" << std::endl; + std::cout << GridLogMessage << "CENSUS 3: power method for sigma_max" << std::endl; + std::cout << GridLogMessage << "=================================================" << std::endl; + MdagMLinearOperator HermOpAdagA(LittleDiracOpPV); + random(CRNG,c_x); + PowerMethod PM; + RealD sigmax2 = PM(HermOpAdagA,c_x); + std::cout << GridLogMessage << "CENSUS: lambda_max(AdagA) = " << sigmax2 + << " sigma_max = " << std::sqrt(sigmax2) << std::endl; + + ////////////////////////////////////////////////////////////////////// + // CENSUS 4: low singular values via Chebyshev-filtered IRL on A^dag A + // The low end of sigma(A_c) is the pseudospectrum of A_c at z=0. + ////////////////////////////////////////////////////////////////////// + std::cout << GridLogMessage << "=================================================" << std::endl; + std::cout << GridLogMessage << "CENSUS 4: Chebyshev-filtered Lanczos, low sigma^2" << std::endl; + std::cout << GridLogMessage << "=================================================" << std::endl; + RealD cheby_hi = (ChebyHi > 0.0) ? ChebyHi : 1.1*sigmax2; + std::cout << GridLogMessage << "Chebyshev filter [" << ChebyLo << "," << cheby_hi << "] order " << ChebyOrder << std::endl; + + // eval/evec/Nconv hoisted out of the block so CENSUS 6 can deflate with them. + std::vector eval(CensusNm); + std::vector evec(CensusNm,Coarse5d); + int Nconv=0; + { + Chebyshev Cheby(ChebyLo,cheby_hi,ChebyOrder); + FunctionHermOp OpCheby(Cheby,HermOpAdagA); + PlainHermOp Op (HermOpAdagA); + + ImplicitlyRestartedLanczos IRL(OpCheby,Op,CensusNstop,CensusNk,CensusNm,CensusTol,CensusMaxIt); + + random(CRNG,c_x); + IRL.calc(eval,evec,c_x,Nconv); + + std::cout << GridLogMessage << "CENSUS: converged " << Nconv << " modes of AdagA" << std::endl; + for(int i=0;i0 ) { +#ifdef HAVE_LIME + std::string evec_file(getenv("CENSUS_EVEC_FILE")); + std::string eval_file = evec_file + ".evals.xml"; + std::cout << GridLogMessage << "CENSUS: saving " << Nconv << " singular vectors to " << evec_file << std::endl; + emptyUserRecord record; + ScidacWriter WR(evec[0].Grid()->IsBoss()); + WR.open(evec_file); + for(int i=0;i eval_out(eval.begin(),eval.begin()+Nconv); // don't shrink shared eval + write(WRx,"evals",eval_out); +#endif + } + } + + // NB: evec/eval stay sized CensusNm (Lattice has no default ctor, so + // std::vector::resize won't instantiate). They match in size, + // which is all DeflatedGuesser asserts; CENSUS 6 only ever indexes [0,Nconv). + + ////////////////////////////////////////////////////////////////////// + // CENSUS 5: half-plane margin from the Hermitian part + // lambda_min(H) = min Re W(A_c) > 0 <=> positive-real (EES applies) + ////////////////////////////////////////////////////////////////////// + std::cout << GridLogMessage << "=================================================" << std::endl; + std::cout << GridLogMessage << "CENSUS 5: Hermitian part H=(A+Adag)/2, half-plane margin" << std::endl; + std::cout << GridLogMessage << "=================================================" << std::endl; + HermitianPartOperator HermPart(LittleDiracOpPV); + + random(CRNG,c_x); + RealD lamHmax = PM(HermPart,c_x); + std::cout << GridLogMessage << "CENSUS: lambda_max(H) = " << lamHmax << std::endl; + + // lambda_min(H): most-negative eigenvalue via Chebyshev-filtered IRL on H. + // A shifted power method cannot separate it from the dense low tail (which is + // why the earlier -0.006 is suspect); Cheby(lo, hi>=lambda_max) amplifies the + // most-negative mode hardest so IRL isolates the true bottom of the spectrum. + RealD hpLo = getenv("HPLANE_CHEBY_LO") ? atof(getenv("HPLANE_CHEBY_LO")) : 0.1; + RealD hpHi = getenv("HPLANE_CHEBY_HI") ? atof(getenv("HPLANE_CHEBY_HI")) : 1.1*lamHmax; + int hpOrder = getenv("HPLANE_CHEBY_ORDER") ? atoi(getenv("HPLANE_CHEBY_ORDER")) : 61; + // Grid's Chebyshev filter MUST be odd order (positive for x < -1, where the low/ + // negative modes map); an even order flips the sign there and the IRL blows up. + if(hpOrder%2==0){ hpOrder++; + std::cout< "< HCheby(hpLo,hpHi,hpOrder); + FunctionHermOp HOpCheby(HCheby,HermPart); + PlainHermOp HOpPlain(HermPart); + ImplicitlyRestartedLanczos HIRL(HOpCheby,HOpPlain,hpNstop,hpNk,hpNm,hpTol,hpMaxIt); + std::vector heval(hpNm); + std::vector hevec(hpNm,Coarse5d); + int hNconv=0; + random(CRNG,c_x); + HIRL.calc(heval,hevec,c_x,hNconv); + RealD lamHmin = (hNconv>0) ? heval[0] : 9.99e99; + for(int kk=0;kk GCR unguaranteed)" << std::endl; + + ////////////////////////////////////////////////////////////////////// + // CENSUS 6: Ritz matrix of the coarse near-null basis + deflated-CG study + // + // C_ij = , S_ij = . + // psi_c are NOT orthonormal (raw near-null projected to coarse), so the + // Rayleigh-Ritz problem is the GENERALISED Hermitian one C v = theta S v. + // Its eigenpairs (theta_i, g_i = sum_j V(j,i) psi_c^j) are the best approximate + // eigenpairs of A^dag A available from span{psi_c}; Eigen normalises so that + // V^dag S V = I, hence = delta_ij and the g_i are an orthonormal + // DeflatedGuesser basis. Compare theta_i to the Lanczos sigma_i^2, then run + // three CG solves on A^dag A: [1] no deflation, [2] Lanczos-eigenvector + // deflated guess, [3] Ritz global-vector deflated guess (g_i treated as pure + // eigenvectors with eigenvalue theta_i). + ////////////////////////////////////////////////////////////////////// + std::cout << GridLogMessage << "=================================================" << std::endl; + std::cout << GridLogMessage << "CENSUS 6: Ritz matrix C_ij = + deflated CG" << std::endl; + std::cout << GridLogMessage << "=================================================" << std::endl; + + std::vector Apsi(nbasis,Coarse5d); + for(int j=0;j(cij.real(),cij.imag()); + Smat(i,j) = std::complex(sij.real(),sij.imag()); + } + } + + { + Eigen::SelfAdjointEigenSolver ses(Smat); + Eigen::MatrixXcd Id = Eigen::MatrixXcd::Identity(nbasis,nbasis); + double SmI = (Smat - Id).norm(); // ||S - I||_F : ~0 iff psi_c orthonormal + std::cout << GridLogMessage << "CENSUS 6: Gram S eig range [" << ses.eigenvalues()(0) + << ", " << ses.eigenvalues()(nbasis-1) + << "] ||S - I||_F = " << SmI + << " (expect ~0: fine vectors are GlobalOrthonormalise'd => psi_c orthonormal)" << std::endl; + } + Eigen::GeneralizedSelfAdjointEigenSolver ges(Cmat,Smat); + Eigen::VectorXd theta = ges.eigenvalues(); // ascending, real + Eigen::MatrixXcd Vr = ges.eigenvectors(); // columns; V^dag S V = I + + int ncmp = std::min((int)nbasis,Nconv); + std::cout << GridLogMessage << "CENSUS 6: Ritz theta vs Lanczos sigma^2 (both ascending):" << std::endl; + for(int i=0;i0) ? std::min(DeflRank,Nconv) : Nconv; + int rankRitz = (DeflRank>0) ? std::min(DeflRank,(int)nbasis) : (int)nbasis; + std::cout << GridLogMessage << "CENSUS 6: CG tol "< CGdefl(CGdeflTol,CGdeflMaxIt,false); + + cg_x = Zero(); + CGdefl(HermOpAdagA,cg_src,cg_x); + std::cout << GridLogMessage << "CENSUS 6: [1] no deflation : iters = " + << CGdefl.IterationsToComplete << " true_resid = " << CGdefl.TrueResidual << std::endl; + + if(rankLanc>0){ + DeflatedGuesser guessL(evec,eval,rankLanc); + guessL(cg_src,cg_x); + CGdefl(HermOpAdagA,cg_src,cg_x); + std::cout << GridLogMessage << "CENSUS 6: [2] Lanczos-evec deflation : iters = " + << CGdefl.IterationsToComplete << " true_resid = " << CGdefl.TrueResidual << std::endl; + } + + { + DeflatedGuesser guessR(gvec,gval,rankRitz); + guessR(cg_src,cg_x); + CGdefl(HermOpAdagA,cg_src,cg_x); + std::cout << GridLogMessage << "CENSUS 6: [3] Ritz-vector deflation : iters = " + << CGdefl.IterationsToComplete << " true_resid = " << CGdefl.TrueResidual << std::endl; + } + + ////////////////////////////////////////////////////////////////////// + // Summary + ////////////////////////////////////////////////////////////////////// + std::cout << GridLogMessage << "=================================================" << std::endl; + std::cout << GridLogMessage << "CENSUS SUMMARY" << std::endl; + std::cout << GridLogMessage << " sigma_max = " << std::sqrt(sigmax2) << std::endl; + std::cout << GridLogMessage << " lambda_max(H) = " << lamHmax << std::endl; + std::cout << GridLogMessage << " lambda_min(H) = " << lamHmin << std::endl; + std::cout << GridLogMessage << " low sigma census : see CENSUS 4 table above" << std::endl; + std::cout << GridLogMessage << " Compare min sigma with |lambda| from Krylov-Schur (Patrick):" << std::endl; + std::cout << GridLogMessage << " sigma_min ~ min|lambda| : effectively normal; deflation rank is the issue" << std::endl; + std::cout << GridLogMessage << " sigma_min << min|lambda|: non-normal; need two-sided/singular-vector deflation" << std::endl; + std::cout << GridLogMessage << "=================================================" << std::endl; + + std::cout << GridLogMessage << "Done" << std::endl; + Grid_finalize(); + return 0; +} diff --git a/examples/Example_pvdagm_halfplane.cc b/examples/Example_pvdagm_halfplane.cc new file mode 100644 index 000000000..f1ce6d28f --- /dev/null +++ b/examples/Example_pvdagm_halfplane.cc @@ -0,0 +1,278 @@ +/* + * Example_pvdagm_halfplane.cc + * + * Standalone fine-operator diagnostic: the EES half-plane margin of the + * (non-Hermitian) PV-preconditioned Mobius DWF operator + * + * A(m_adj) = D_adj^dag D_light (D_adj plays the Pauli-Villars role) + * + * as a function of the adjoint mass m_adj, dialled from the light quark mass + * up to the Pauli-Villars mass (=1). No coarse grid, no subspace, no Lanczos + * -- pure power-method spectral tests on the fine grid. + * + * Purpose: A is the LEFT preconditioner for inverting the light operator. + * To solve D_light X = B we iterate the preconditioned system + * (D_adj^dag D_light) X = D_adj^dag B , + * whose solution X is independent of m_adj -- only the conditioning and the + * iterative convergence change. m_adj = m_light is the usual CGNR (symmetric + * normal equations); m_adj = 1 is the Pauli-Villars preconditioned system. + * The sweep asks which m_adj keeps the preconditioned operator well-behaved + * (positive-real / EES-guaranteed) while buying the wider spectral range. + * + * For the Hermitian part H(A) = (A + A^dag)/2 we measure, per m_adj: + * + * lambda_max(H) -- power method on H + * lambda_min(H) -- power method on (sI - H) => min Re W(A), the half-plane + * margin. EES (Eisenstat-Elman-Schultz 1983, Thm 3.3) + * GUARANTEES GCR convergence with rate + * [ 1 - lambda_min(H)^2 / sigma_max^2 ]^{1/2} + * ONLY when lambda_min(H) > 0 (positive-real / A's field + * of values in the open right half-plane). A negative + * value means the guarantee is lost (not that GCR + * diverges); the magnitude is then the distance-to- + * positive-realness, i.e. the shift/deflation needed to + * recover it. + * sigma_max -- power method on A^dag A (= A.HermOp) + * + * Endpoints: + * m_adj = m_light => A = M^dag M, Hermitian PD, positive-real by + * construction, lambda_min(H) = sigma_min^2 > 0 (the + * squared / CGNR operator). + * m_adj = 1 => A = PV^dag M, the standard PVdagM operator. + * + * Env: MASS, M5, MOBIUS_B, MOBIUS_C, LS, CONFIG, + * MADJ_LIST (comma separated) OR MADJ_MIN / MADJ_MAX / MADJ_N (geometric). + * + * Caveat: lambda_min(H) via a shifted power method can be soft when it sits + * near zero over a dense low spectrum. The SIGN and the TREND across m_adj + * are the robust signal; confirm an individual near-zero value with a proper + * shifted Lanczos if it is load-bearing. + */ + +#include + +using namespace std; +using namespace Grid; + +////////////////////////////////////////////////////////////////////// +// A = PV^dag M : Op = _PV.Mdag . _Mat.M , AdjOp = _Mat.Mdag . _PV.M +////////////////////////////////////////////////////////////////////// +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){ // A^dag A + Field tmp(in.Grid()); + Op(in,tmp); + AdjOp(tmp,out); + } +}; + +////////////////////////////////////////////////////////////////////// +// H = (A + A^dag)/2 for a general non-Hermitian LinearOperator A. +////////////////////////////////////////////////////////////////////// +template +class HermitianPartLinOp : public LinearOperatorBase { + LinearOperatorBase &_A; +public: + HermitianPartLinOp(LinearOperatorBase &A): _A(A) {}; + 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){ HermOp(in,out); } + void AdjOp (const Field &in, Field &out){ HermOp(in,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()); + _A.Op(in,out); // A in + _A.AdjOp(in,tmp); // A^dag in + out = 0.5*(out + tmp); + } +}; + +////////////////////////////////////////////////////////////////////// +// s*I - Op : power method on this gives s - lambda_min(Op) for Hermitian Op. +////////////////////////////////////////////////////////////////////// +template +class ShiftedNegatedOperator : public LinearOperatorBase { + LinearOperatorBase &_Op; + RealD s; +public: + ShiftedNegatedOperator(RealD _s, LinearOperatorBase &Op): _Op(Op), s(_s) {}; + 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){ HermOp(in,out); } + void AdjOp (const Field &in, Field &out){ HermOp(in,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){ + _Op.HermOp(in,out); + out = s*in - out; + } +}; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + RealD mass = 0.00078; + RealD M5 = 1.8; + RealD b = 1.5; + RealD c = 0.5; + int Ls = 24; + std::string config("ckpoint_lat.1000"); + + if(getenv("MASS")) mass = atof(getenv("MASS")); + if(getenv("M5")) M5 = atof(getenv("M5")); + if(getenv("MOBIUS_B")) b = atof(getenv("MOBIUS_B")); + if(getenv("MOBIUS_C")) c = atof(getenv("MOBIUS_C")); + if(getenv("LS")) Ls = atoi(getenv("LS")); + if(getenv("CONFIG")) config = std::string(getenv("CONFIG")); + + // Adjoint-mass sweep: explicit list, or geometric MADJ_MIN..MADJ_MAX in MADJ_N steps. + std::vector madj_list; + if(getenv("MADJ_LIST")){ + std::stringstream ss(getenv("MADJ_LIST")); + std::string tok; + while(std::getline(ss,tok,',')) if(tok.size()) madj_list.push_back(std::stod(tok)); + } else { + int N = getenv("MADJ_N") ? atoi(getenv("MADJ_N")) : 6; + RealD lo = getenv("MADJ_MIN") ? atof(getenv("MADJ_MIN")) : mass; + RealD hi = getenv("MADJ_MAX") ? atof(getenv("MADJ_MAX")) : 1.0; + GRID_ASSERT(N>=1); + for(int i=0;i auto + int HalfChebyOrder = getenv("HALF_CHEBY_ORDER") ? atoi(getenv("HALF_CHEBY_ORDER")) : 61; + // Grid's Chebyshev filter MUST be odd order: only then is the polynomial positive + // for x < -1, the region the low/negative modes map to. An even order flips the + // sign there, the filtered operator explodes negative, and the IRL never converges. + if(HalfChebyOrder%2==0){ HalfChebyOrder++; + std::cout< "< lat = {48,48,48,96}; + + GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(lat, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); + GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); + GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); + + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers({5,6,7,8}); + + std::cout << GridLogMessage << "PARAM: MASS(light) " << mass << " M5 " << M5 + << " b " << b << " c " << c << " Ls " << Ls << std::endl; + std::cout << GridLogMessage << "PARAM: CONFIG " << config << std::endl; + + LatticeGaugeField Umu(UGrid); + FieldMetaData header; + std::cout << GridLogMessage << "Reading gauge field " << config << std::endl; + NerscIO::readConfiguration(Umu,header,config); + + // Fixed light operator (never changes across the sweep). + MobiusFermionD Dlight(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid, mass, M5, b, c); + + LatticeFermionD x(FGrid); + + std::cout << GridLogMessage << "=================================================" << std::endl; + std::cout << GridLogMessage << "FINE HALF-PLANE SWEEP A(m_adj) = D_adj^dag D_light" << std::endl; + std::cout << GridLogMessage << " m_adj = " << mass << " => M^dag M (positive-real); m_adj = 1 => PVdagM" << std::endl; + std::cout << GridLogMessage << "=================================================" << std::endl; + + for(auto madj : madj_list){ + + MobiusFermionD Dadj(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid, madj, M5, b, c); + + PVdagMLinearOperator A(Dlight,Dadj); // A = Dadj^dag Dlight + HermitianPartLinOp H(A); + + PowerMethod PM; + + random(RNG5,x); RealD lamHmax = PM(H,x); + + // lambda_min(H): most-negative eigenvalue via Chebyshev-filtered IRL on H. + // Cheby(lo,hi) amplifies eigenvalues below lo; with hi>=lambda_max(H) the most + // negative mode is amplified hardest, so IRL isolates the true bottom of the + // (possibly indefinite) spectrum where the shifted power method could not. + RealD fhi = (HalfChebyHi>0.0)? HalfChebyHi : 1.1*lamHmax; + Chebyshev Cheby(HalfChebyLo,fhi,HalfChebyOrder); + FunctionHermOp OpCheby(Cheby,H); + PlainHermOp OpPlain(H); + ImplicitlyRestartedLanczos IRL(OpCheby,OpPlain,HalfNstop,HalfNk,HalfNm,HalfTol,HalfMaxIt); + std::vector heval(HalfNm); + std::vector hevec(HalfNm,FGrid); + int hNconv=0; + random(RNG5,x); + IRL.calc(heval,hevec,x,hNconv); + RealD lamHmin = (hNconv>0) ? heval[0] : 9.99e99; + for(int kk=0;kk 0.0); + RealD ratefac = posreal ? std::sqrt(1.0 - lamHmin*lamHmin/sigmax2) : 0.0; // EES per-iter + RealD iters8 = (posreal && ratefac < 1.0) ? std::log(1.0e-8)/std::log(ratefac) : 0.0; + + std::cout << GridLogMessage << "HALFPLANE: m_adj " << madj + << " lambda_min(H) " << lamHmin + << " lambda_max(H) " << lamHmax + << " sigma_max " << sigmax + << " positive_real " << (posreal ? "YES" : "NO ") + << (posreal + ? (" EES_rate " + std::to_string(ratefac) + " EES_iters(1e-8) " + std::to_string(iters8)) + : (" margin_below_zero " + std::to_string(-lamHmin) + " (EES guarantee lost)")) + << std::endl; + } + + std::cout << GridLogMessage << "=================================================" << std::endl; + std::cout << GridLogMessage << "Reading: lambda_min(H) > 0 => EES guarantees GCR at the quoted rate." << std::endl; + std::cout << GridLogMessage << " crossing to < 0 as m_adj -> 1 marks loss of positive-realness." << std::endl; + std::cout << GridLogMessage << " (non-normality: eigenvalues may still be right-half-plane.)" << std::endl; + std::cout << GridLogMessage << "Done" << std::endl; + + Grid_finalize(); + return 0; +} diff --git a/examples/Example_pvdagm_mrhs_3level.cc b/examples/Example_pvdagm_mrhs_3level.cc new file mode 100644 index 000000000..582dfb138 --- /dev/null +++ b/examples/Example_pvdagm_mrhs_3level.cc @@ -0,0 +1,613 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./examples/Example_pvdagm_mrhs_3level.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. + + See the full license in the file "LICENSE" in the top level distribution directory + *************************************************************************************/ + /* END LEGAL */ + +// MultiRHS (valence) THREE-level multigrid for PVdagM. +// +// This is exactly the plain three-level algorithm of Example_pvdagm_3level_SVDdefl.cc +// with L3_DEFL=0 (NO deflation), applied to the enlarged block-diagonal mRHS system: +// the coarse and coarse-coarse levels run a SINGLE Krylov (one GCR polynomial, inner +// products summed over rhs) on the packed 6D mrhs fields, so both coarse levels batch +// through GEMM (MultiGeneralCoarsenedMatrix) -- the valence throughput win at BOTH levels. +// +// Level structure (each coarse level is a single-field PGCR on a packed 6D mrhs field): +// L1 (fine) : std::vector, MrhsPGCRNonHermitian on PVdagM, +// preconditioned by the L1->L2 mrhs V-cycle (MrhsTwoLevelMG). +// L2 (coarse) : 6D mrhs coarse field, PGCR, preconditioned by the L2->L3 mrhs +// V-cycle (MrhsCoarseThreeLevelPrec) -- coarse-coarse correction + coarse smoother. +// L3 (coarse-coarse): 6D mrhs coarse-coarse field, PGCR (the innermost solve). +// +// RAW-NULL DISCIPLINE (critical -- see project_block_orthogonalise_leak): the L2->L3 +// aggregation MUST be built from RAW fine near-null vectors (pre block-GS). We take a +// raw copy of the loaded subspace BEFORE the L1->L2 CoarsenOperator (which block- +// orthonormalises in place) and project THAT. Guards print || - I||: ~0.23 = +// content preserved, ~N_coarse = the e_k leak is back. +// +// Env: MASS SUBSPACE_FILE NRHS +// BLOCK (dotted, default 2.2.2.2) BLOCK2 (dotted, default 2.2.3.3) +// FineSmootherShift FineSmootherOrder +// CoarseSmootherShift CoarseSmootherNstep +// CoarseSolverTol CoarseSolverOrder +// L3_TOL L3_MAXIT L3_NSTEP +// OuterMmax OuterNstep OuterTol + +#include +#include +#include +#include + +using namespace std; +using namespace Grid; + +RealD FineSmootherShift = 0.1; +int FineSmootherOrder = 16; +RealD CoarseSmootherShift = 0.1; +int CoarseSmootherNstep = 4; +RealD CoarseSolverTol = 0.03; +int CoarseSolverOrder = 200; +RealD L3Tol = 2.5e-1; +int L3MaxIt = 50; +int L3Nstep = 50; +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("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("L3_TOL")) L3Tol = atof(getenv("L3_TOL")); + if(getenv("L3_MAXIT")) L3MaxIt = atoi(getenv("L3_MAXIT")); + if(getenv("L3_NSTEP")) L3Nstep = atoi(getenv("L3_NSTEP")); + 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: CoarseSmootherShift" << CoarseSmootherShift<< std::endl; + std::cout << GridLogMessage << "PARAM: CoarseSmootherNstep" << CoarseSmootherNstep<< std::endl; + std::cout << GridLogMessage << "PARAM: CoarseSolverTol " << CoarseSolverTol << std::endl; + std::cout << GridLogMessage << "PARAM: CoarseSolverOrder " << CoarseSolverOrder << std::endl; + std::cout << GridLogMessage << "PARAM: L3_TOL " << L3Tol << 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 + 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), and shifted variant for smoothers. +////////////////////////////////////////////////////////////////////// +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); } +}; + +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); } +}; + +// Generic shift wrapper (for the coarse-level smoother on the 6D mrhs coarse operator). +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); } +}; + +////////////////////////////////////////////////////////////////////// +// mrhs interfaces + single-polynomial mrhs PGCR (verbatim from Example_pvdagm_mrhs.cc): +// reductions summed over rhs -> one alpha/beta per step for the enlarged system. +////////////////////////////////////////////////////////////////////// +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; // 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; } + 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); + 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); return cp; + } +}; + +////////////////////////////////////////////////////////////////////// +// L2->L3 mrhs V-cycle: a LinearFunction on the 6D mrhs COARSE field. +// Mirrors Example_pvdagm_mrhs.cc's MrhsTwoLevelMG one level down, and the +// single-RHS MGPreconditioner of Example_pvdagm_3level_SVDdefl.cc: +// out = in (trivial pre) +// r = in - A_coarse out +// restrict (unpack 6D coarse -> blockProject -> pack 6D coarse-coarse) +// ONE coarse-coarse solve (L3, GEMM) +// prolong (unpack -> blockPromote -> pack); out += correction +// r = in - A_coarse out +// coarse smoother (shifted 6D coarse op); out += smooth(r) +////////////////////////////////////////////////////////////////////// +template +class MrhsCoarseThreeLevelPrec : public LinearFunction { +public: + LinearOperatorBase &_CoarseOp; // mrhs coarse op (6D) + LinearFunction &_CoarseSmoother; // shifted 6D coarse smoother + MultiRHSBlockProject &_Projector; // L2->L3 (vector-based) + LinearFunction &_CoarseCoarseSolve; // L3 solve (6D cc) + 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; double t; + CoarseField vec1(in.Grid()); + CoarseField vec2(in.Grid()); + + // trivial pre-smoother + out = in; + + // residual (6D coarse) + _CoarseOp.Op(out,vec1); sub(vec1,in,vec1); + + // restrict: unpack 6D coarse -> vector -> blockProject -> vector -> pack 6D cc + std::vector csplit(nrhs,_Coarse5d); + std::vector ccsplit(nrhs,_CoarseCoarse5d); + CoarseCoarseField CCsrc(_CoarseCoarseMrhs); + CoarseCoarseField CCsol(_CoarseCoarseMrhs); + + t=-usecond(); + for(int r=0;rL3 restrict took "< blockPromote -> pack 6D coarse; add correction + t=-usecond(); + for(int r=0;rL3 prolong took "<L2 mrhs V-cycle (verbatim from Example_pvdagm_mrhs.cc): +// per-rhs fine smoother + batched restriction + ONE coarse solve + batched prolong. +// The coarse solve passed in is now itself three-level (preconditioned by L2->L3). +////////////////////////////////////////////////////////////////////// +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(); double t; + std::vector vec1(nrhs,fgrid),vec2(nrhs,fgrid); + for(int r=0;r Csrc_split(nrhs,_CoarseGrid), Csol_split(nrhs,_CoarseGrid); + CoarseVector CsrcMrhs(_CoarseGridMrhs), CsolMrhs(_CoarseGridMrhs); + t=-usecond(); + _Projector.blockProject(vec1,Csrc_split); + for(int r=0;r 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); + + // 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; + + // Level 2 blocking (default 2,2,3,3) -- matches Example_pvdagm_3level_SVDdefl + Coordinate cclatt = clatt; + Coordinate Block2({2,2,3,3}); + if ( getenv("BLOCK2") ){ GridCmdOptionIntVector(std::string(getenv("BLOCK2")),Block2); GRID_ASSERT(Block2.size()==4); } + for(int d=0;d<4;d++){ GRID_ASSERT(clatt[d]%Block2[d]==0); cclatt[d]=clatt[d]/Block2[d]; } + std::cout << GridLogMessage << "Block2 " << Block2 << " coarse-coarse lattice " << cclatt << std::endl; + + GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); + GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d); + GridCartesian *CoarseCoarse4d = SpaceTimeGrid::makeFourDimGrid(cclatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); + GridCartesian *CoarseCoarse5d = SpaceTimeGrid::makeFiveDimGrid(1,CoarseCoarse4d); + + // 6D mrhs grids: rhs is dim 0, SIMD across rhs (pattern: Test_general_coarse_hdcg_phys48.cc) + Coordinate mpi=GridDefaultMpi(); + Coordinate rhMpi ({1,1,mpi[0],mpi[1],mpi[2],mpi[3]}); + Coordinate rhSimd({vComplex::Nsimd(),1,1,1,1,1}); + Coordinate rhLatt ({nrhs,1,clatt[0], clatt[1], clatt[2], clatt[3]}); + Coordinate rhLatt2({nrhs,1,cclatt[0],cclatt[1],cclatt[2],cclatt[3]}); + GridCartesian *CoarseMrhs = new GridCartesian(rhLatt, rhSimd,rhMpi); + GridCartesian *CoarseCoarseMrhs = new GridCartesian(rhLatt2,rhSimd,rhMpi); + + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers({5,6,7,8}); + + LatticeGaugeField Umu(UGrid); + std::cout << GridLogMessage << "Reading gauge field" << std::endl; + FieldMetaData header; + std::string file("/ccs/home/poare/ckpoint_lat.1000"); + 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; + typedef ShiftedPVdagMLinearOperator ShiftedPVdagM_t; + + // Level 1 tensor types + typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef MultiGeneralCoarsenedMatrix MrhsLittleDiracOperator; + typedef LittleDiracOperator::CoarseVector CoarseVector; + typedef Aggregation Subspace; + + // Level 2 tensor types (coarsening deepens the nest by one iScalar -- see CLAUDE.md) + typedef CoarseVector::vector_object CoarseSiteObj; + typedef iScalar vTTComplex; + typedef GeneralCoarsenedMatrix LittleDiracOperatorL2; + typedef MultiGeneralCoarsenedMatrix MrhsLittleDiracOperatorL2; + typedef LittleDiracOperatorL2::CoarseVector CoarseCoarseVector; + typedef Aggregation SubspaceL2; + + PVdagM_t PVdagM(Ddwf,Dpv); + ShiftedPVdagM_t ShiftedPVdagM(FineSmootherShift,Ddwf,Dpv); + + NextToNearestStencilGeometry5D geom (Coarse5d); + NextToNearestStencilGeometry5D geom2(CoarseCoarse5d); // 33-point at L2->L3, matching SVDdefl + + ////////////////////////////////////////////////////////////////////// + // Subspace: load RAW (no Orthogonalise!), 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 (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 of the fine null vectors BEFORE CoarsenOperator block-orthonormalises in place. + std::vector rawNull(nbasis,FGrid); + for(int k=0;kL2 and L2->L3 with SINGLE-RHS machinery, import into the mrhs + // operators via CopyMatrix. The single-RHS L1->L2 coarse operator must stay + // alive to be the "fine" operator for the L2->L3 coarsening, so BOTH single-RHS + // ops (and their padded _A) live in one scope and free together. [MEMORY: this + // is the setup peak -- L1->L2 padded _A (~large at 2^4) + L2->L3 padded _A.] + ////////////////////////////////////////////////////////////////////// + MrhsLittleDiracOperator mrhsLittleDiracOpPV(geom, CoarseMrhs); + MrhsLittleDiracOperatorL2 mrhsLittleDiracOpL2(geom2, CoarseCoarseMrhs); + MultiRHSBlockProject MrhsProjector; + MultiRHSBlockProject MrhsProjectorL2; + { + // --- L1->L2 single-RHS coarse operator (kept alive for the L2->L3 coarsening) --- + LittleDiracOperator LittleDiracOpPV(geom,FGrid,Coarse5d); + LittleDiracOpPV.CoarsenOperator(PVdagM, AggregatesGCR); // orthonormalises AggregatesGCR.subspace in place + mrhsLittleDiracOpPV.CopyMatrix(LittleDiracOpPV); + MrhsProjector.Allocate(nbasis,FGrid,Coarse5d); + MrhsProjector.ImportBasis(AggregatesGCR.subspace); // orthonormalised, matches the coarse op + NonHermitianLinearOperator LinOpCoarse(LittleDiracOpPV); + + // --- psi_coarse = P^dag (RAW fine null) -> Galerkin images, NOT e_k --- + std::vector psi_coarse(nbasis,Coarse5d); + for(int k=0;k - I||_F = "<L3 single-RHS coarsening (coarsen the single-RHS LinOpCoarse) --- + SubspaceL2 AggregatesL2(CoarseCoarse5d,Coarse5d,cb); + for(int k=0;k psi_cc(nbasis,CoarseCoarse5d); + for(int k=0;k - I||_F = "< mrhsLinOpCoarse(mrhsLittleDiracOpPV); + NonHermitianLinearOperator mrhsLinOpCC(mrhsLittleDiracOpL2); + + ////////////////////////////////////////////////////////////////////// + // Solvers, innermost first. + ////////////////////////////////////////////////////////////////////// + TrivialPrecon simpleC; + TrivialPrecon simpleCC; + TrivialPrecon simple_fine; + + // L3 (coarse-coarse) solve: PGCR on the 6D cc operator + PrecGeneralisedConjugateResidualNonHermitian + L3PGCR(L3Tol,L3MaxIt,mrhsLinOpCC,simpleCC,L3Nstep,L3Nstep); + L3PGCR.Level(3); + L3PGCR.Name("CCouter"); + L3PGCR.SetZeroGuess(1); // caller zeroes CCsol + + // L2 coarse smoother: shifted 6D coarse op, fixed nstep + ShiftedLinearOperator ShiftedMrhsCoarse(CoarseSmootherShift, mrhsLinOpCoarse); + PrecGeneralisedConjugateResidualNonHermitian + CoarseSmootherGCR(0.01,1,ShiftedMrhsCoarse,simpleC,CoarseSmootherNstep,CoarseSmootherNstep); + CoarseSmootherGCR.Level(2); + CoarseSmootherGCR.Name("Csmoother"); + CoarseSmootherGCR.SetZeroGuess(1); // caller zeroes vec2 + + // L2->L3 V-cycle preconditioner (operates on 6D coarse field) + MrhsCoarseThreeLevelPrec + L2to3Precon(mrhsLinOpCoarse, CoarseSmootherGCR, MrhsProjectorL2, L3PGCR, + Coarse5d, CoarseCoarse5d, CoarseCoarseMrhs, nrhs); + + // L2 coarse solve: PGCR on 6D coarse op, preconditioned by the L2->L3 V-cycle + PrecGeneralisedConjugateResidualNonHermitian + L2PGCR(CoarseSolverTol, CoarseSolverOrder/16, mrhsLinOpCoarse, L2to3Precon, 16, 16); + L2PGCR.Level(2); + L2PGCR.Name("Couter"); + L2PGCR.SetZeroGuess(1); // caller zeroes CsolMrhs + + // Fine smoother (per-rhs, looped in the L1->L2 V-cycle) + PrecGeneralisedConjugateResidualNonHermitian + SmootherGCR(0.0,1,ShiftedPVdagM,simple_fine,FineSmootherOrder,FineSmootherOrder); + SmootherGCR.Level(1); + SmootherGCR.Name("Fsmoother"); + SmootherGCR.SetZeroGuess(1); // caller zeroes vec2[r] + + // L1->L2 V-cycle (fine); its coarse solve is the three-level L2PGCR + typedef PrecGeneralisedConjugateResidualNonHermitian FineSmoother_t; + MrhsTwoLevelMG + ThreeLevelPrecon(PVdagM, SmootherGCR, MrhsProjector, L2PGCR, Coarse5d, CoarseMrhs); + + // Outer mrhs solve + MrhsPGCRNonHermitian + L1PGCR(OuterTol,1000,PVdagM,ThreeLevelPrecon,OuterMmax,OuterNstep); + L1PGCR.Level(1); + L1PGCR.Name("Fouter"); + L1PGCR.SetZeroGuess(1); // sol[r]=Zero() at source setup + + ////////////////////////////////////////////////////////////////////// + // Sources and solve + ////////////////////////////////////////////////////////////////////// + std::vector src(nrhs,FGrid), 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. + + See the full license in the file "LICENSE" in the top level distribution directory + *************************************************************************************/ + /* END LEGAL */ + +// MultiRHS (valence) THREE-level multigrid for PVdagM with a DENSE, EXACT, +// non-iterative coarse-coarse bottom == Example_pvdagm_mrhs_3level.cc with +// the L3 PGCR replaced by DistributedDenseInverse (the machinery validated in +// Example_pvdagm_3level_dense.cc: VERIFY ~6e-5, apply ~6ms single-RHS). +// +// The two headline results COMPOSED: mrhs GEMM batching at the coarse level +// (the ~12x valence win) x supercoarse-dense exact bottom (the coarsest level +// reduced from dominant cost to noise). Under mrhs the dense apply becomes a +// GEMM against the resident row-slab: ONE GlobalSum of nrhs packed vectors + +// ONE slab pass for all rhs -- per-rhs bottom cost BELOW the single-RHS 6ms. +// +// Level structure: +// L1 (fine) : std::vector, MrhsPGCRNonHermitian on PVdagM, +// preconditioned by the L1->L2 mrhs V-cycle (MrhsTwoLevelMG). +// L2 (coarse) : 6D mrhs coarse field, PGCR, preconditioned by the L2->L3 mrhs +// V-cycle -- DENSE coarse-coarse correction + coarse smoother. +// L3 (coarse-coarse): DENSE row-distributed A^{-1} (cgetrf_64 + blocked cgetrs_64, +// fp32, ILP64 -- rank N = gSites(CC) x nbasis, e.g. 69120 at +// BLOCK2=8.4.3.6 x nb60). DENSE_CC=0 reverts to the L3 PGCR. +// +// SUPERCOARSE default: BLOCK2 = 8,4,3,6 (CC = [3,6,8,8], the dense floor on +// --mpi 3.6.4.4; sigma spectrum is FLAT across all 60 vectors => keep nb60). +// +// RAW-NULL DISCIPLINE unchanged (see project_block_orthogonalise_leak): guards +// print || - I||_F at both levels (~0.23 good, ~sqrt(N) = e_k leak). +// +// Build (Frontier/HIP): LIBS += -lrocsolver -lrocblas. +// Run: --device-mem 40000 (solve-phase residency; EvictAll makes the setup +// window; boss GCD holds 38.2GB dA + 141MB dB at BLOCK2=8.4.3.6 nb60). +// +// Env: MASS SUBSPACE_FILE NRHS +// BLOCK (dotted, default 2.2.2.2) BLOCK2 (dotted, default 8.4.3.6) +// FineSmootherShift FineSmootherOrder +// CoarseSmootherShift CoarseSmootherNstep +// CoarseSolverTol CoarseSolverOrder +// DENSE_CC (default 1) DENSE_CC_CHECK +// L3_TOL L3_MAXIT L3_NSTEP (iterative branch only) +// OuterMmax OuterNstep OuterTol + +#include +#include +#include +#include + +#include +#include + +#ifdef GRID_HIP +#include +#endif + +using namespace std; +using namespace Grid; + +RealD FineSmootherShift = 0.1; +int FineSmootherOrder = 16; +RealD CoarseSmootherShift = 0.1; +int CoarseSmootherNstep = 4; +RealD CoarseSolverTol = 0.03; +int CoarseSolverOrder = 200; +RealD L3Tol = 2.5e-1; +int L3MaxIt = 50; +int L3Nstep = 50; +RealD OuterTol = 1.0e-8; +int OuterMmax = 8; +int OuterNstep = 8; +int Nrhs = 12; +int UseDenseCC = 1; +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("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("L3_TOL")) L3Tol = atof(getenv("L3_TOL")); + if(getenv("L3_MAXIT")) L3MaxIt = atoi(getenv("L3_MAXIT")); + if(getenv("L3_NSTEP")) L3Nstep = atoi(getenv("L3_NSTEP")); + 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")); + if(getenv("DENSE_CC")) UseDenseCC = atoi(getenv("DENSE_CC")); + + std::cout << GridLogMessage << "PARAM: MASS " << mass << std::endl; + std::cout << GridLogMessage << "PARAM: NRHS " << Nrhs << std::endl; + std::cout << GridLogMessage << "PARAM: DENSE_CC " << UseDenseCC << std::endl; + std::cout << GridLogMessage << "PARAM: FineSmootherShift " << FineSmootherShift << std::endl; + std::cout << GridLogMessage << "PARAM: FineSmootherOrder " << FineSmootherOrder << std::endl; + std::cout << GridLogMessage << "PARAM: CoarseSmootherShift" << CoarseSmootherShift<< std::endl; + std::cout << GridLogMessage << "PARAM: CoarseSmootherNstep" << CoarseSmootherNstep<< std::endl; + std::cout << GridLogMessage << "PARAM: CoarseSolverTol " << CoarseSolverTol << std::endl; + std::cout << GridLogMessage << "PARAM: CoarseSolverOrder " << CoarseSolverOrder << std::endl; + std::cout << GridLogMessage << "PARAM: L3_TOL " << L3Tol << 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 + 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), and shifted variant for smoothers. +////////////////////////////////////////////////////////////////////// +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); } +}; + +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); } +}; + +// Generic shift wrapper (for the coarse-level smoother on the 6D mrhs coarse operator). +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); } +}; + +////////////////////////////////////////////////////////////////////////////////////// +// DistributedDenseInverse: exact solve of a (small) coarse operator by explicit, +// row-distributed dense inverse. VALIDATED single-RHS in Example_pvdagm_3level_dense.cc +// (VERIFY ~6e-5). Setup: probe columns -> exact chunked GlobalSum gather streamed +// into the boss-GCD device buffer -> cgetrf_64 (ILP64) -> rows of A^{-1} produced +// blockwise via cgetrs_64 on identity-column blocks and broadcast; each rank keeps +// the rows of ITS OWN sites (ownership-aligned => apply needs ONE GlobalSum only). +// +// NEW HERE: ApplyBatch -- the mrhs form. ONE GlobalSum of all nrhs packed vectors, +// ONE pass over the slab computing all rhs (GEMM shape). Per-rhs cost < single-RHS. +// +// Tensor-depth agnostic: site scalar_object treated as nbasis contiguous ComplexD. +////////////////////////////////////////////////////////////////////////////////////// +template +class DistributedDenseInverse : public LinearFunction { +public: + using LinearFunction::operator(); + typedef typename Field::vector_object vobj; + typedef typename vobj::scalar_object sobj; + + GridBase *grid; + LinearOperatorBase &_Op; // kept for env-gated defect checks (DENSE_CC_CHECK) + int nbasis; + int64_t N; // dense rank = gSites * nbasis + int nd; + int lsites; // my local sites + int nrows; // my rows = lsites * nbasis + std::vector myLcoor; // local coordinate of my site ss + std::vector myGsite; // global lex site index of my site ss + std::vector slab; // nrows x N row-major rows of A^{-1} + // CHUNKROWS=1024: fatter cgetrs_64 calls (trsm efficiency grows with nrhs; + // 256 gave ~7s/chunk => 32min setup) and 4x fewer broadcasts. Buffers 566MB. + static const int64_t CHUNKROWS = 1024; + static const int MRHS_MAX = 32; +#ifdef GRID_HIP + rocblas_float_complex *dSlab = nullptr; // device-resident copy of the slab (all ranks) + rocblas_handle applyHandle; // per-rank handle for the apply GEMM + // Persistent apply buffers, hoisted out of the per-call path: pinned host + // staging (full PCIe rate, no page faults) + device X/Y (no per-call + // hipMalloc/hipFree). Sized once for MRHS_MAX. + ComplexF *hXapply = nullptr; + ComplexF *hYapply = nullptr; + rocblas_float_complex *dXapply = nullptr; + rocblas_float_complex *dYapply = nullptr; + int devSum = 0; // DENSE_DEVICE_SUM=1: allreduce the DEVICE buffer + // (GPU-aware MPI, no host staging); default host sum +#endif + // NB: the apply GEMM (Y = slab^dag X) is a tiny-output/huge-K shape that + // under-fills the GPU (~13ms). The fix is a software split-K via + // GridBLAS.gemmBatched (see MultiRHSBlockCGLinalg.h / 2409.03904 Fig 11) — + // NOT a raw rocblas strided-batched batch, which hung on Frontier and was + // removed. TODO: reimplement through GridBLAS when the ~1.4 s/RHS is wanted. + + DistributedDenseInverse(LinearOperatorBase &Op, GridBase *g, int nbasis_) + : grid(g), _Op(Op), nbasis(nbasis_) + { + GRID_ASSERT( sizeof(sobj) == nbasis*sizeof(ComplexD) ); // site object == nbasis ComplexD + nd = grid->_ndimension; + N = grid->gSites() * nbasis; + lsites = grid->lSites(); + nrows = lsites * nbasis; + + Coordinate ldims = grid->LocalDimensions(); + Coordinate gdims = grid->GlobalDimensions(); + + std::cout << GridLogMessage << "DistributedDenseInverse: N = " << N + << " (" << grid->gSites() << " sites x " << nbasis << ")" + << " rows/rank = " << nrows + << " slab = " << (double)nrows*N*sizeof(ComplexF)/1024./1024. << " MB/rank" + << std::endl; + + //////////////////////////////////////////////////////////////////// + // Enumerate my sites: local coords and global lexicographic indices + //////////////////////////////////////////////////////////////////// + myLcoor.resize(lsites); + myGsite.resize(lsites); + for(int ss=0; ss_lstart[d] + lcoor[d]; + int64_t gsite; + Lexicographic::IndexFromCoor(gcoor, gsite, gdims); + myLcoor[ss] = lcoor; + myGsite[ss] = gsite; + } + + slab.resize((uint64_t)nrows * N); + + double t0 = usecond(); + //////////////////////////////////////////////////////////////////// + // 0. Slab cache: SLAB_FILE= -> per-rank raw file .. + // Present: load, skip probe/gather/factor/solve (setup ~free on + // resubmits and sweeps). Absent: full setup, then write it. + // VERIFY below runs in BOTH paths, certifying loaded data. + //////////////////////////////////////////////////////////////////// + bool loaded = false; + char *sfile = getenv("SLAB_FILE"); + std::string slabfile; + if (sfile) { + slabfile = std::string(sfile) + "." + std::to_string(grid->ThisRank()); + FILE *f = fopen(slabfile.c_str(),"rb"); + if (f) { + int64_t hdr[4] = {0,0,0,0}; + GRID_ASSERT( fread(hdr,sizeof(int64_t),4,f) == 4 ); + GRID_ASSERT( hdr[0] == (int64_t)0x44454E5345 ); // magic "DENSE" + GRID_ASSERT( hdr[1] == N && hdr[2] == (int64_t)nrows && hdr[3] == (int64_t)nbasis ); + uint64_t nelem = (uint64_t)nrows * N; + GRID_ASSERT( fread(&slab[0], sizeof(ComplexF), nelem, f) == nelem ); + fclose(f); + loaded = true; + std::cout << GridLogMessage << "DistributedDenseInverse: slab loaded from " + << slabfile << " -- skipping probe/factor/solve" << std::endl; + } else { + std::cout << GridLogMessage << "DistributedDenseInverse: slab cache " << slabfile + << " absent -- full setup, will write it" << std::endl; + } + } + if (!loaded) { + //////////////////////////////////////////////////////////////////// + // 1. PROBE assembly of my rows of A: column (jsite,b) = Op(e_{jsite,b}) + //////////////////////////////////////////////////////////////////// + Field e(grid); + Field Ae(grid); + Coordinate gcoorj(nd); + int64_t gsitesN = grid->gSites(); + for(int64_t jsite=0; jsiteIsBoss(); + std::vector Afull; +#ifdef GRID_HIP + rocblas_float_complex *dA = nullptr; + rocblas_float_complex *dB = nullptr; // getrs_64 RHS block (CHUNKROWS identity columns) + rocblas_handle rochandle; + int64_t *dIpiv = nullptr; // ILP64 pivots, live from factor to last getrs + uint64_t Abytes = (uint64_t)N * N * sizeof(ComplexF); + // Make HBM space for the naked hipMalloc: flush the device-copy layer. + // (FreePool of the allocator free-list awaits the type-dispatched fix.) + MemoryManager::EvictAll(); + // MemoryManager::FreePool(); + if (boss) { + auto aerr = hipMalloc((void **)&dA, Abytes); + if (aerr != hipSuccess) { + std::cout << GridLogMessage << "DistributedDenseInverse: hipMalloc of " + << Abytes/1024./1024./1024. << " GB FAILED -- reduce --device-mem " + << "or enable the fixed MemoryManager::FreePool()" << std::endl; + GRID_ASSERT(aerr == hipSuccess); + } + std::cout << GridLogMessage << "DistributedDenseInverse: device inversion buffer allocated (" + << Abytes/1024./1024./1024. << " GB)" << std::endl; + } +#else + if (boss) Afull.resize((uint64_t)N * N); +#endif + { + std::unordered_map rowmap; // global row -> my slab row + for(int ss=0; ss chunk((uint64_t)CHUNKROWS * N); + for(int64_t row0=0; row0second) * N; + uint64_t dst = (uint64_t)(r-row0) * N; + for(int64_t j=0;jGlobalSumVector(&chunk[0], (int)nelem); + if (boss) { +#ifdef GRID_HIP + GRID_ASSERT( hipMemcpy((char *)dA + (uint64_t)row0*N*sizeof(ComplexF), + &chunk[0], nelem*sizeof(ComplexF), + hipMemcpyHostToDevice) == hipSuccess ); +#else + uint64_t dst = (uint64_t)row0 * N; + for(uint64_t i=0;i LU of A^T. + //////////////////////////////////////////////////////////////////// + if (boss) { +#ifdef GRID_HIP + std::cout << GridLogMessage << "DistributedDenseInverse: rocSOLVER cgetrf_64 (ILP64 LU) N=" << N + << " in place on resident device buffer" << std::endl; + auto hst = rocblas_create_handle(&rochandle); + std::cout << GridLogMessage << "DistributedDenseInverse: rocblas handle status " << (int)hst << std::endl; + int64_t *dInfo; + GRID_ASSERT( hipMalloc((void **)&dIpiv, N*sizeof(int64_t)) == hipSuccess ); + GRID_ASSERT( hipMalloc((void **)&dInfo, sizeof(int64_t)) == hipSuccess ); + auto st1 = rocsolver_cgetrf_64(rochandle, (int64_t)N, (int64_t)N, dA, (int64_t)N, dIpiv, dInfo); + hipDeviceSynchronize(); + int64_t info_h = -1; + hipMemcpy(&info_h, dInfo, sizeof(int64_t), hipMemcpyDeviceToHost); + std::cout << GridLogMessage << "DistributedDenseInverse: cgetrf_64 status " << (int)st1 + << " info = " << (int)info_h << std::endl; + GRID_ASSERT(st1 == rocblas_status_success); + GRID_ASSERT(info_h == 0); + hipFree(dInfo); + GRID_ASSERT( hipMalloc((void **)&dB, (uint64_t)CHUNKROWS*N*sizeof(ComplexF)) == hipSuccess ); + // dA holds the LU of A^T; rows of A^{-1} are produced blockwise below via + // cgetrs_64 on identity-column blocks: A^T X = E => X columns = rows of + // A^{-1}, in exactly the linear layout the harvest expects. +#else + // Eigen fallback: small local CPU tests only. + std::cout << GridLogMessage << "DistributedDenseInverse: Eigen fallback inversion N=" << N + << (N > 10000 ? " (WARNING: SLOW; use the HIP/rocSOLVER path)" : "") + << std::endl; + typedef Eigen::Matrix,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> MatF; + Eigen::Map A(reinterpret_cast*>(&Afull[0]), N, N); + MatF Ainv = A.inverse(); + A = Ainv; +#endif + } + double t3 = usecond(); + std::cout << GridLogMessage << "DistributedDenseInverse: factorisation took " + << (t3-t2)/1.0e6 << " s" << std::endl; + + //////////////////////////////////////////////////////////////////// + // 4. Blocked solve + broadcast: rows of A^{-1} chunk by chunk; each + // rank keeps the rows of its own sites (ownership-aligned). + //////////////////////////////////////////////////////////////////// + { + std::unordered_map rowmap; + for(int ss=0; ss chunk((uint64_t)CHUNKROWS * N); + for(int64_t row0=0; row0Broadcast(0, &chunk[0], nelem*sizeof(ComplexF)); + for(int64_t r=row0; rsecond) * N; + uint64_t src = (uint64_t)(r-row0) * N; + for(int64_t j=0;j x((uint64_t)N, ComplexD(0.0,0.0)); + for(int ss=0; ssGlobalSumVector(&x[0], (int)N); + + std::vector y(nrows); + thread_for(r, nrows, { + ComplexD acc(0.0,0.0); + const ComplexF *row = &slab[(uint64_t)r * N]; + for(int64_t j=0; j X(nX, ComplexD(0.0,0.0)); + { + autoView(iv, in, CpuRead); + Coordinate c6(nd+1); + for(int ss=0; ssGlobalSumVector(&X[0], (int)nX); + std::vector Y((uint64_t)nrows*nr); + thread_for(r_, nrows, { + ComplexD acc[MRHS_MAX]; + for(int rr=0; rr +class MrhsDenseCCSolve : public LinearFunction { +public: + DistributedDenseInverse &_Dense; + GridBase *_CoarseCoarse5d; + int _nrhs; + MrhsDenseCCSolve(DistributedDenseInverse &D, GridBase *cc5d, int nrhs) + : _Dense(D), _CoarseCoarse5d(cc5d), _nrhs(nrhs) {} + using LinearFunction::operator(); + virtual void operator()(const CoarseCoarseField &in, CoarseCoarseField &out){ + if ( getenv("DENSE_CC_CHECK") ) { + // Audit path: per-rhs 5D unpack so ApplyBatch can run the _Op defect + // check per rhs. ~50ms/call of slice/split overhead -- audit only. + CoarseCoarseField tmp(in.Grid()); + tmp = in; + std::vector split_in (_nrhs,_CoarseCoarse5d); + std::vector split_out(_nrhs,_CoarseCoarse5d); + for(int r=0;r<_nrhs;r++) ExtractSliceFast(split_in[r], tmp, r, 0); + _Dense.ApplyBatch(split_in, split_out); + for(int r=0;r<_nrhs;r++) InsertSliceFast(split_out[r], out, r, 0); + } else { + _Dense.ApplyBatch6D(in, out, _nrhs); + } + } +}; + +////////////////////////////////////////////////////////////////////// +// mrhs interfaces + single-polynomial mrhs PGCR (verbatim from Example_pvdagm_mrhs.cc) +////////////////////////////////////////////////////////////////////// +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; // 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; } + 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); + 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); return cp; + } +}; + +////////////////////////////////////////////////////////////////////// +// L2->L3 mrhs V-cycle: LinearFunction on the 6D mrhs COARSE field. +// The coarse-coarse solve slot now takes EITHER the dense mrhs solve +// (DENSE_CC=1) or the L3 PGCR (DENSE_CC=0). +////////////////////////////////////////////////////////////////////// +template +class MrhsCoarseThreeLevelPrec : public LinearFunction { +public: + LinearOperatorBase &_CoarseOp; // mrhs coarse op (6D) + LinearFunction &_CoarseSmoother; // shifted 6D coarse smoother + MultiRHSBlockProject &_Projector; // L2->L3 (vector-based) + LinearFunction &_CoarseCoarseSolve; // L3 solve (6D cc) + 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; double t; + CoarseField vec1(in.Grid()); + CoarseField vec2(in.Grid()); + + // trivial pre-smoother + out = in; + + // residual (6D coarse) + _CoarseOp.Op(out,vec1); sub(vec1,in,vec1); + + // restrict: unpack 6D coarse -> vector -> blockProject -> vector -> pack 6D cc + std::vector csplit(nrhs,_Coarse5d); + std::vector ccsplit(nrhs,_CoarseCoarse5d); + CoarseCoarseField CCsrc(_CoarseCoarseMrhs); + CoarseCoarseField CCsol(_CoarseCoarseMrhs); + + t=-usecond(); + for(int r=0;rL3 restrict took "< blockPromote -> pack 6D coarse; add correction + t=-usecond(); + for(int r=0;rL3 prolong took "<L2 mrhs V-cycle (verbatim from Example_pvdagm_mrhs.cc) +////////////////////////////////////////////////////////////////////// +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(); double t; + std::vector vec1(nrhs,fgrid),vec2(nrhs,fgrid); + for(int r=0;r Csrc_split(nrhs,_CoarseGrid), Csol_split(nrhs,_CoarseGrid); + CoarseVector CsrcMrhs(_CoarseGridMrhs), CsolMrhs(_CoarseGridMrhs); + t=-usecond(); + _Projector.blockProject(vec1,Csrc_split); + for(int r=0;r 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); + + // 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; + + // Level 2 blocking: SUPERCOARSE default 8,4,3,6 -> CC [3,6,8,8], the dense floor. + Coordinate cclatt = clatt; + Coordinate Block2({8,4,3,6}); + if ( getenv("BLOCK2") ){ GridCmdOptionIntVector(std::string(getenv("BLOCK2")),Block2); GRID_ASSERT(Block2.size()==4); } + for(int d=0;d<4;d++){ GRID_ASSERT(clatt[d]%Block2[d]==0); cclatt[d]=clatt[d]/Block2[d]; } + std::cout << GridLogMessage << "Block2 " << Block2 << " coarse-coarse lattice " << cclatt << std::endl; + + GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); + GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d); + GridCartesian *CoarseCoarse4d = SpaceTimeGrid::makeFourDimGrid(cclatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); + GridCartesian *CoarseCoarse5d = SpaceTimeGrid::makeFiveDimGrid(1,CoarseCoarse4d); + + // 6D mrhs grids: rhs is dim 0, SIMD across rhs + Coordinate mpi=GridDefaultMpi(); + Coordinate rhMpi ({1,1,mpi[0],mpi[1],mpi[2],mpi[3]}); + Coordinate rhSimd({vComplex::Nsimd(),1,1,1,1,1}); + Coordinate rhLatt ({nrhs,1,clatt[0], clatt[1], clatt[2], clatt[3]}); + Coordinate rhLatt2({nrhs,1,cclatt[0],cclatt[1],cclatt[2],cclatt[3]}); + GridCartesian *CoarseMrhs = new GridCartesian(rhLatt, rhSimd,rhMpi); + GridCartesian *CoarseCoarseMrhs = new GridCartesian(rhLatt2,rhSimd,rhMpi); + + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers({5,6,7,8}); + + LatticeGaugeField Umu(UGrid); + std::cout << GridLogMessage << "Reading gauge field" << std::endl; + FieldMetaData header; + std::string file("/ccs/home/poare/ckpoint_lat.1000"); + 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; + typedef ShiftedPVdagMLinearOperator ShiftedPVdagM_t; + + // Level 1 tensor types + typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef MultiGeneralCoarsenedMatrix MrhsLittleDiracOperator; + typedef LittleDiracOperator::CoarseVector CoarseVector; + typedef Aggregation Subspace; + + // Level 2 tensor types (coarsening deepens the nest by one iScalar) + typedef CoarseVector::vector_object CoarseSiteObj; + typedef iScalar vTTComplex; + typedef GeneralCoarsenedMatrix LittleDiracOperatorL2; + typedef MultiGeneralCoarsenedMatrix MrhsLittleDiracOperatorL2; + typedef LittleDiracOperatorL2::CoarseVector CoarseCoarseVector; + typedef Aggregation SubspaceL2; + + PVdagM_t PVdagM(Ddwf,Dpv); + ShiftedPVdagM_t ShiftedPVdagM(FineSmootherShift,Ddwf,Dpv); + + NextToNearestStencilGeometry5D geom (Coarse5d); + NextToNearestStencilGeometry5D geom2(CoarseCoarse5d); + + ////////////////////////////////////////////////////////////////////// + // Subspace: load RAW (no Orthogonalise!), 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 (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 of the fine null vectors BEFORE CoarsenOperator block-orthonormalises in place. + std::vector rawNull(nbasis,FGrid); + for(int k=0;kL2 and L2->L3 with SINGLE-RHS machinery; import to mrhs via + // CopyMatrix. The L2 (coarse-coarse) single-RHS operator is HOISTED to + // main scope: it is tiny at the supercoarse blocking and the dense setup + // (probe/VERIFY/defect checks) needs it alive for the whole run. + ////////////////////////////////////////////////////////////////////// + MrhsLittleDiracOperator mrhsLittleDiracOpPV(geom, CoarseMrhs); + MrhsLittleDiracOperatorL2 mrhsLittleDiracOpL2(geom2, CoarseCoarseMrhs); + MultiRHSBlockProject MrhsProjector; + MultiRHSBlockProject MrhsProjectorL2; + + LittleDiracOperatorL2 LittleDiracOpL2(geom2,Coarse5d,CoarseCoarse5d); + NonHermitianLinearOperator LinOpCC5d(LittleDiracOpL2); + + { + // --- L1->L2 single-RHS coarse operator (scoped: its padded _A is the memory peak) --- + LittleDiracOperator LittleDiracOpPV(geom,FGrid,Coarse5d); + LittleDiracOpPV.CoarsenOperator(PVdagM, AggregatesGCR); // orthonormalises AggregatesGCR.subspace in place + mrhsLittleDiracOpPV.CopyMatrix(LittleDiracOpPV); + MrhsProjector.Allocate(nbasis,FGrid,Coarse5d); + MrhsProjector.ImportBasis(AggregatesGCR.subspace); + NonHermitianLinearOperator LinOpCoarse(LittleDiracOpPV); + + // --- psi_coarse = P^dag (RAW fine null) -> Galerkin images, NOT e_k --- + std::vector psi_coarse(nbasis,Coarse5d); + for(int k=0;k - I||_F = "<L3 single-RHS coarsening --- + SubspaceL2 AggregatesL2(CoarseCoarse5d,Coarse5d,cb); + for(int k=0;k psi_cc(nbasis,CoarseCoarse5d); + for(int k=0;k - I||_F = "< mrhsLinOpCoarse(mrhsLittleDiracOpPV); + NonHermitianLinearOperator mrhsLinOpCC(mrhsLittleDiracOpL2); + + ////////////////////////////////////////////////////////////////////// + // DENSE coarse-coarse bottom (constructed AFTER the fine coarsening frees + // its memory peak; needs only the hoisted single-RHS LinOpCC5d). + ////////////////////////////////////////////////////////////////////// + std::unique_ptr> DenseCC; + std::unique_ptr> MrhsDenseCC; + if (UseDenseCC) { + std::cout << GridLogMessage << "**********************************************" << std::endl; + std::cout << GridLogMessage << " Dense CC inverse setup (mrhs bottom)" << std::endl; + std::cout << GridLogMessage << "**********************************************" << std::endl; + DenseCC.reset(new DistributedDenseInverse(LinOpCC5d, CoarseCoarse5d, nbasis)); + MrhsDenseCC.reset(new MrhsDenseCCSolve(*DenseCC, CoarseCoarse5d, nrhs)); + } + + ////////////////////////////////////////////////////////////////////// + // Solvers, innermost first. + ////////////////////////////////////////////////////////////////////// + TrivialPrecon simpleC; + TrivialPrecon simpleCC; + TrivialPrecon simple_fine; + + // L3 (coarse-coarse) iterative solve: PGCR on the 6D cc operator (DENSE_CC=0 branch) + PrecGeneralisedConjugateResidualNonHermitian + L3PGCR(L3Tol,L3MaxIt,mrhsLinOpCC,simpleCC,L3Nstep,L3Nstep); + L3PGCR.Level(3); + L3PGCR.Name("CCouter"); + + LinearFunction *ccSolve; + if (UseDenseCC) ccSolve = MrhsDenseCC.get(); + else ccSolve = &L3PGCR; + + // L2 coarse smoother: shifted 6D coarse op, fixed nstep + ShiftedLinearOperator ShiftedMrhsCoarse(CoarseSmootherShift, mrhsLinOpCoarse); + PrecGeneralisedConjugateResidualNonHermitian + CoarseSmootherGCR(0.01,1,ShiftedMrhsCoarse,simpleC,CoarseSmootherNstep,CoarseSmootherNstep); + CoarseSmootherGCR.Level(2); + CoarseSmootherGCR.Name("Csmoother"); + CoarseSmootherGCR.SetZeroGuess(1); // caller zeroes vec2: skip r0 apply every L2 iteration + + // L2->L3 V-cycle preconditioner (operates on 6D coarse field) + MrhsCoarseThreeLevelPrec + L2to3Precon(mrhsLinOpCoarse, CoarseSmootherGCR, MrhsProjectorL2, *ccSolve, + Coarse5d, CoarseCoarse5d, CoarseCoarseMrhs, nrhs); + + // L2 coarse solve: PGCR on 6D coarse op, preconditioned by the L2->L3 V-cycle + PrecGeneralisedConjugateResidualNonHermitian + L2PGCR(CoarseSolverTol, CoarseSolverOrder/16, mrhsLinOpCoarse, L2to3Precon, 16, 16); + L2PGCR.Level(2); + L2PGCR.Name("Couter"); + L2PGCR.SetZeroGuess(1); // caller zeroes CsolMrhs; restarts still recompute r + + // Fine smoother (per-rhs, looped in the L1->L2 V-cycle) + PrecGeneralisedConjugateResidualNonHermitian + SmootherGCR(0.0,1,ShiftedPVdagM,simple_fine,FineSmootherOrder,FineSmootherOrder); + SmootherGCR.Level(1); + SmootherGCR.Name("Fsmoother"); + SmootherGCR.SetZeroGuess(1); // caller zeroes vec2[r]: saves 12 fine mults/outer + + // L1->L2 V-cycle (fine); its coarse solve is the three-level L2PGCR + typedef PrecGeneralisedConjugateResidualNonHermitian FineSmoother_t; + MrhsTwoLevelMG + ThreeLevelPrecon(PVdagM, SmootherGCR, MrhsProjector, L2PGCR, Coarse5d, CoarseMrhs); + + // Outer mrhs solve + MrhsPGCRNonHermitian + L1PGCR(OuterTol,1000,PVdagM,ThreeLevelPrecon,OuterMmax,OuterNstep); + L1PGCR.Level(1); + L1PGCR.Name("Fouter"); + L1PGCR.SetZeroGuess(1); // sol[r]=Zero() below; restarts recompute r as always + + ////////////////////////////////////////////////////////////////////// + // Sources and solve + ////////////////////////////////////////////////////////////////////// + std::vector src(nrhs,FGrid), sol(nrhs,FGrid); + for(int r=0;r