portelli/Grid
1
0 Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2026-10-08 08:48:06 +01:00
Code Releases Activity

Synch multigrid rewrite to repository and

Browse Source 
broad unicode elimination effort
This commit is contained in:
Peter Boyle committed 2026-09-29 12:36:34 -04:00
1 parent 865c338dd2
commit 3bd5e8883f
73 files changed
+1018 -3219
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

examples/Example_mdagm.cc
+2 -2
View File
@@ -87,7 +87,7 @@ int main (int argc, char ** argv)
HermOpAdaptor<LatticeFermionD> HermFineOp(MdagMOp);
// ── Coarse geometry ────────────────────────────────────────────────────
typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef DeprecatedGeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef LittleDiracOperator::CoarseVector CoarseVector;
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
@@ -171,7 +171,7 @@ int main (int argc, char ** argv)
LatticeFermionD src(FGrid); random(RNG5, src);
LatticeFermionD result(FGrid); result = Zero();
TwoLevelADEF2<LatticeFermionD, CoarseVector, Subspace>
DeprecatedTwoLevelADEF2<LatticeFermionD, CoarseVector, Subspace>
HDCG(1.0e-8, 1000,
HermFineOp,
Smoother,
examples/Example_mdagm_cg.cc
+2 -2
View File
@@ -31,7 +31,7 @@ Author: Peter Boyle <pboyle@bnl.gov>
//
// Operator hierarchy:
// Fine: M†M, acted on via HermOpAdaptor so Op = HermOp = M†M
// Coarse: 33-point GeneralCoarsenedMatrix (NextToNearest, 2-hop M†M)
// Coarse: 33-point DeprecatedGeneralCoarsenedMatrix (NextToNearest, 2-hop M†M)
//
// Setup:
// 1. Chebyshev filter (T_600 x T_2500) on M†M to build near-null subspace {ψᵢ}
@@ -198,7 +198,7 @@ int main(int argc, char **argv)
///////////////////////////////////////////////////////////
// Coarse operator: 33-point stencil for 2-hop M†M
///////////////////////////////////////////////////////////
typedef GeneralCoarsenedMatrix<vSpinColourVector, vTComplex, nbasis> LittleDiracOperator;
typedef DeprecatedGeneralCoarsenedMatrix<vSpinColourVector, vTComplex, nbasis> LittleDiracOperator;
typedef LittleDiracOperator::CoarseVector CoarseVector;
NextToNearestStencilGeometry5D geom(Coarse5d);
examples/Example_pvdagm.cc
+2 -2
View File
@@ -364,7 +364,7 @@ void runMG(
) {
// typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
// typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
// typedef DeprecatedGeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
// typedef LittleDiracOperator::CoarseVector CoarseVector;
ParseEnvironment();
@@ -620,7 +620,7 @@ int main (int argc, char ** argv)
// assert(nbasis <= Nevecs); // need to have enough evecs
typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef DeprecatedGeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef LittleDiracOperator::CoarseVector CoarseVector;
NextToNearestStencilGeometry5D geom(Coarse5d);
examples/Example_pvdagm_3level.cc
+23 -23
View File
@@ -190,8 +190,8 @@ public:
}
};
// Lüscher deflated guesser (arXiv:0706.2298 Sec A.3) for a non-Hermitian solve.
// C_{st} = <psi[s] | LinOp | psi[t]>; guess = sum_s c_s psi[s] where c = C^{-1} psi† src.
// Luscher deflated guesser (arXiv:0706.2298 Sec A.3) for a non-Hermitian solve.
// C_{st} = <psi[s] | LinOp | psi[t]>; guess = sum_s c_s psi[s] where c = C^{-1} psi^dag src.
template<class Field>
class LuscherGuesser : public LinearFunction<Field> {
const std::vector<Field> &psi;
@@ -342,7 +342,7 @@ void runMG(
TrivialPrecon<LatticeFermionD> simple_fine;
//////////////////////////////////////////////////////////////////////
// Level 0→1: coarsen PVdagM, build LinOpCoarse
// Level 0->1: coarsen PVdagM, build LinOpCoarse
//////////////////////////////////////////////////////////////////////
LittleDiracOperator LittleDiracOpPV(geom, FGrid, Coarse5d);
LittleDiracOpPV.CoarsenOperator(PVdagM, AggregatesPD);
@@ -366,7 +366,7 @@ void runMG(
//////////////////////////////////////////////////////////////////////
// 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.
// Used as the aggregation basis for Level 1->2 coarsening.
//////////////////////////////////////////////////////////////////////
std::vector<CoarseVector> psi_coarse(nbasis, Coarse5d);
for (int k = 0; k < nbasis; k++)
@@ -396,22 +396,22 @@ void runMG(
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 << "Coarse null matrix ||C - C^dag||/||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.
// 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).
// geometry 3x6x4x4 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<vTComplex>, so CComplex
// for the L1→L2 level must be iScalar<vTComplex>, not vTComplex.
// for the L1->L2 level must be iScalar<vTComplex>, not vTComplex.
typedef typename CoarseVector::vector_object CoarseSiteObj;
typedef iScalar<vTComplex> vTTComplex;
typedef GeneralCoarsenedMatrix<CoarseSiteObj,vTTComplex,NB> LittleDiracOperatorL2;
typedef DeprecatedGeneralCoarsenedMatrix<CoarseSiteObj,vTTComplex,NB> LittleDiracOperatorL2;
typedef typename LittleDiracOperatorL2::CoarseVector CoarseCoarseVector;
typedef Aggregation<CoarseSiteObj,vTTComplex,NB> SubspaceL2;
typedef MGPreconditioner<CoarseSiteObj,vTTComplex,NB> L1to2MG;
@@ -429,12 +429,12 @@ void runMG(
TrivialPrecon<CoarseCoarseVector> simpleCC;
//////////////////////////////////////////////////////////////////////
// Lüscher deflation guesser for L3PGCR.
// Luscher deflation guesser for L3PGCR.
// Step 1: project psi_coarse[k] (promoted fine null vectors) to
// CoarseCoarseVector space — these cover the zero-momentum
// 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
// 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} = <psi_cc[s]|LinOpCC|psi_cc[t]> over the
// full augmented basis and invert directly via Eigen LU.
@@ -476,7 +476,7 @@ void runMG(
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;
std::cout << GridLogMessage << "Coarse-coarse deflation matrix ||Ccc-Ccc^dag||/||Ccc|| = " << normCccmCdag/normCcc << std::endl;
}
Eigen::MatrixXcd Ccc_inv = Ccc.inverse();
LuscherGuesser<CoarseCoarseVector> CCDeflGuesser(psi_cc, Ccc_inv);
@@ -489,7 +489,7 @@ void runMG(
L3PGCR.Name("CCouter");
//////////////////////////////////////////////////////////////////////
// Coarse-level GCR smoother for Level 1→2 V-cycle.
// 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.
//////////////////////////////////////////////////////////////////////
@@ -505,7 +505,7 @@ void runMG(
CoarseSmootherGCR.Name("Csmoother");
//////////////////////////////////////////////////////////////////////
// Level 1→2 V-cycle preconditioner.
// Level 1->2 V-cycle preconditioner.
//////////////////////////////////////////////////////////////////////
L1to2MG L1to2Precon(AggregatesL2,
LinOpCoarse,
@@ -513,7 +513,7 @@ void runMG(
CoarseSmootherGCR, // post-smoother: 12 GCR steps
LinOpCC,
L3PGCR,
CCDeflGuesser); // Lüscher guesser: psi_cc C^{-1} psi_cc†
CCDeflGuesser); // Luscher guesser: psi_cc C^{-1} psi_cc^dag
//////////////////////////////////////////////////////////////////////
// Standalone Level 1 two-level solve test.
@@ -548,7 +548,7 @@ void runMG(
f_src = one;
// Pre-smoother: none (TrivialPrecon); post-smoother: shifted PGCR.
// Coarse solver: L2MGsolver (PGCR preconditioned by Level 1→2 V-cycle).
// Coarse solver: L2MGsolver (PGCR preconditioned by Level 1->2 V-cycle).
TwoLevelMG ThreeLevelPrecon(AggregatesPD,
PVdagM,
simple_fine,
@@ -592,7 +592,7 @@ int main (int argc, char ** argv)
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)
// Level 1 coarse grid: block 2^4 from fine (48x48x48x96 -> 24x24x24x48, 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;
@@ -600,13 +600,13 @@ int main (int argc, char ** argv)
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).
// Level 2 coarse-coarse grid: block 2,2,3,3 from Level 1 (24x24x24x48 -> 12x12x8x16, 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.
// 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). ✓
// 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;
@@ -635,7 +635,7 @@ int main (int argc, char ** argv)
typedef PVdagMLinearOperator<MobiusFermionD,LatticeFermionD> PVdagM_t;
typedef ShiftedPVdagMLinearOperator<MobiusFermionD,LatticeFermionD> ShiftedPVdagM_t;
typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef DeprecatedGeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef LittleDiracOperator::CoarseVector CoarseVector;
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
typedef MGPreconditioner<vSpinColourVector,vTComplex,nbasis> TwoLevelMG;
examples/Example_pvdagm_3level_dense.cc
+12 -12
View File
@@ -212,8 +212,8 @@ public:
}
};
// Lüscher deflated guesser (arXiv:0706.2298 Sec A.3) for a non-Hermitian solve.
// C_{st} = <psi[s] | LinOp | psi[t]>; guess = sum_s c_s psi[s] where c = C^{-1} psi† src.
// Luscher deflated guesser (arXiv:0706.2298 Sec A.3) for a non-Hermitian solve.
// C_{st} = <psi[s] | LinOp | psi[t]>; guess = sum_s c_s psi[s] where c = C^{-1} psi^dag src.
template<class Field>
class LuscherGuesser : public LinearFunction<Field> {
const std::vector<Field> &psi;
@@ -704,7 +704,7 @@ void runMG(
TrivialPrecon<LatticeFermionD> simple_fine;
//////////////////////////////////////////////////////////////////////
// Level 0→1: coarsen PVdagM, build LinOpCoarse
// Level 0->1: coarsen PVdagM, build LinOpCoarse
//////////////////////////////////////////////////////////////////////
LittleDiracOperator LittleDiracOpPV(geom, FGrid, Coarse5d);
LittleDiracOpPV.CoarsenOperator(PVdagM, AggregatesPD);
@@ -736,7 +736,7 @@ void runMG(
// Optional sigma-ordering of psi_coarse (SVD_REORDER set): replace the crude
// first-NB_CC slice with the NB_CC most-null directions of span(psi_coarse)
// under LinOpCoarse. Nullness measure = singular values (eig of the
// Gram-whitened Psi†A†APsi), NOT the numerical range Q†AQ which
// Gram-whitened Psi^dag A^dag APsi), NOT the numerical range Q^dag AQ which
// non-normality contaminates. The printed sigma spectrum shows where the
// truncation cliff sits. Unset => raw first-30 (crude GS-ordered slice).
//////////////////////////////////////////////////////////////////////
@@ -787,7 +787,7 @@ void runMG(
}
//////////////////////////////////////////////////////////////////////
// Level 1→2: SUPERCOARSE aggregation using psi_coarse as subspace.
// Level 1->2: SUPERCOARSE aggregation using psi_coarse as subspace.
// Maximal block {8,4,3,6}: CC = [3,6,8,8] = the dense-invertible floor.
//
// UNBLOCKING TRUNCATION NB_CC = 30 (first-30 slice of psi_coarse):
@@ -811,7 +811,7 @@ void runMG(
typedef typename CoarseVector::vector_object CoarseSiteObj;
typedef iScalar<vTComplex> vTTComplex;
typedef GeneralCoarsenedMatrix<CoarseSiteObj,vTTComplex,NB_CC> LittleDiracOperatorL2;
typedef DeprecatedGeneralCoarsenedMatrix<CoarseSiteObj,vTTComplex,NB_CC> LittleDiracOperatorL2;
typedef typename LittleDiracOperatorL2::CoarseVector CoarseCoarseVector;
typedef Aggregation<CoarseSiteObj,vTTComplex,NB_CC> SubspaceL2;
typedef MGPreconditioner<CoarseSiteObj,vTTComplex,NB_CC> L1to2MG;
@@ -830,7 +830,7 @@ void runMG(
//////////////////////////////////////////////////////////////////////
// CC solve: DENSE (exact, non-iterative) by default; DENSE_CC=0 gives
// the previous iterative L3PGCR + Lüscher-guesser path for A/B.
// the previous iterative L3PGCR + Luscher-guesser path for A/B.
//////////////////////////////////////////////////////////////////////
int use_dense = 1;
if (getenv("DENSE_CC")) use_dense = atoi(getenv("DENSE_CC"));
@@ -859,7 +859,7 @@ void runMG(
ccGuess = &simpleCC; // exact solve ignores/overwrites any guess
} else {
////////////////////////////////////////////////////////////////////
// Lüscher deflation guesser (arXiv:0706.2298 A.3) for the iterative CC
// Luscher deflation guesser (arXiv:0706.2298 A.3) for the iterative CC
// solve, as in the earlier supercoarse configuration.
////////////////////////////////////////////////////////////////////
psi_cc.resize(nbasis, CoarseCoarse5d);
@@ -900,7 +900,7 @@ void runMG(
}
//////////////////////////////////////////////////////////////////////
// Coarse-level GCR smoother for Level 1→2 V-cycle.
// Coarse-level GCR smoother for Level 1->2 V-cycle.
//////////////////////////////////////////////////////////////////////
RealD coarse_smoother_shift = 0.1;
int coarse_smoother_nstep = 2;
@@ -914,7 +914,7 @@ void runMG(
CoarseSmootherGCR.SetZeroGuess(1); // post-smoother slot: caller zeroes vec2 (NOT L2MGsolver: it takes the Luscher guess)
//////////////////////////////////////////////////////////////////////
// Level 1→2 V-cycle preconditioner.
// Level 1->2 V-cycle preconditioner.
//////////////////////////////////////////////////////////////////////
L1to2MG L1to2Precon(AggregatesL2,
LinOpCoarse,
@@ -998,7 +998,7 @@ int main (int argc, char ** argv)
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)
// Level 1 coarse grid: block 2^4 from fine (48x48x48x96 -> 24x24x24x48, Ls=1)
Coordinate clatt = lat_size;
// Coordinate Block1({2,2,2,2});
// Coordinate Block2({8,4,3,6});
@@ -1047,7 +1047,7 @@ int main (int argc, char ** argv)
typedef PVdagMLinearOperator<MobiusFermionD,LatticeFermionD> PVdagM_t;
typedef ShiftedPVdagMLinearOperator<MobiusFermionD,LatticeFermionD> ShiftedPVdagM_t;
typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef DeprecatedGeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef LittleDiracOperator::CoarseVector CoarseVector;
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
typedef MGPreconditioner<vSpinColourVector,vTComplex,nbasis> TwoLevelMG;
examples/Example_pvdagm_3level_madj.cc
+17 -17
View File
@@ -323,7 +323,7 @@ void runMG(
TrivialPrecon<LatticeFermionD> simple_fine;
//////////////////////////////////////////////////////////////////////
// Level 0→1: coarsen PVdagM, build LinOpCoarse
// Level 0->1: coarsen PVdagM, build LinOpCoarse
//////////////////////////////////////////////////////////////////////
LittleDiracOperator LittleDiracOpPV(geom, FGrid, Coarse5d);
LittleDiracOpPV.CoarsenOperator(PVdagM, AggregatesPD);
@@ -348,7 +348,7 @@ void runMG(
//////////////////////////////////////////////////////////////////////
// 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.
// Used as the aggregation basis for Level 1->2 coarsening.
//////////////////////////////////////////////////////////////////////
std::vector<CoarseVector> psi_coarse(nbasis, Coarse5d);
for (int k = 0; k < nbasis; k++)
@@ -378,22 +378,22 @@ void runMG(
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 << "Coarse null matrix ||C - C^dag||/||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.
// 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).
// geometry 3x6x4x4 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<vTComplex>, so CComplex
// for the L1→L2 level must be iScalar<vTComplex>, not vTComplex.
// for the L1->L2 level must be iScalar<vTComplex>, not vTComplex.
typedef typename CoarseVector::vector_object CoarseSiteObj;
typedef iScalar<vTComplex> vTTComplex;
typedef GeneralCoarsenedMatrix<CoarseSiteObj,vTTComplex,NB> LittleDiracOperatorL2;
typedef DeprecatedGeneralCoarsenedMatrix<CoarseSiteObj,vTTComplex,NB> LittleDiracOperatorL2;
typedef typename LittleDiracOperatorL2::CoarseVector CoarseCoarseVector;
typedef Aggregation<CoarseSiteObj,vTTComplex,NB> SubspaceL2;
typedef MGPreconditioner<CoarseSiteObj,vTTComplex,NB> L1to2MG;
@@ -412,14 +412,14 @@ void runMG(
// Level 2 solver: plain GCR, no further coarsening
//////////////////////////////////////////////////////////////////////
TrivialPrecon<CoarseCoarseVector> simpleCC;
// L3PGCR is an inner solver inside the L1→2 V-cycle; does not need to converge
// 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<CoarseCoarseVector> L3PGCR(1.0e-4,5,LinOpCC,simpleCC,64,64);
L3PGCR.Level(3);
//////////////////////////////////////////////////////////////////////
// Coarse-level GCR smoother for Level 1→2 V-cycle.
// 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.
//////////////////////////////////////////////////////////////////////
@@ -433,7 +433,7 @@ void runMG(
CoarseSmootherGCR.Level(2);
//////////////////////////////////////////////////////////////////////
// Level 1→2 V-cycle preconditioner.
// Level 1->2 V-cycle preconditioner.
//////////////////////////////////////////////////////////////////////
L1to2MG L1to2Precon(AggregatesL2,
LinOpCoarse,
@@ -474,7 +474,7 @@ void runMG(
f_src = one;
// Pre-smoother: none (TrivialPrecon); post-smoother: shifted PGCR.
// Coarse solver: L2MGsolver (PGCR preconditioned by Level 1→2 V-cycle).
// Coarse solver: L2MGsolver (PGCR preconditioned by Level 1->2 V-cycle).
TwoLevelMG ThreeLevelPrecon(AggregatesPD,
PVdagM,
simple_fine,
@@ -518,7 +518,7 @@ int main (int argc, char ** argv)
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)
// Level 1 coarse grid: block 2^4 from fine (48x48x48x96 -> 24x24x24x48, 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;
@@ -526,13 +526,13 @@ int main (int argc, char ** argv)
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).
// Level 2 coarse-coarse grid: block 2,2,3,3 from Level 1 (24x24x24x48 -> 12x12x8x16, 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.
// 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). ✓
// 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;
@@ -566,7 +566,7 @@ int main (int argc, char ** argv)
typedef PVdagMLinearOperator<MobiusFermionD,LatticeFermionD> PVdagM_t;
typedef ShiftedPVdagMLinearOperator<MobiusFermionD,LatticeFermionD> ShiftedPVdagM_t;
typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef DeprecatedGeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef LittleDiracOperator::CoarseVector CoarseVector;
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
typedef MGPreconditioner<vSpinColourVector,vTComplex,nbasis> TwoLevelMG;
examples/Example_pvdagm_4level.cc
+29 -29
View File
@@ -190,8 +190,8 @@ public:
}
};
// Lüscher deflated guesser (arXiv:0706.2298 Sec A.3) for a non-Hermitian solve.
// C_{st} = <psi[s] | LinOp | psi[t]>; guess = sum_s c_s psi[s] where c = C^{-1} psi† src.
// Luscher deflated guesser (arXiv:0706.2298 Sec A.3) for a non-Hermitian solve.
// C_{st} = <psi[s] | LinOp | psi[t]>; guess = sum_s c_s psi[s] where c = C^{-1} psi^dag src.
template<class Field>
class LuscherGuesser : public LinearFunction<Field> {
const std::vector<Field> &psi;
@@ -343,7 +343,7 @@ void runMG(
TrivialPrecon<LatticeFermionD> simple_fine;
//////////////////////////////////////////////////////////////////////
// Level 0→1: coarsen PVdagM, build LinOpCoarse
// Level 0->1: coarsen PVdagM, build LinOpCoarse
//////////////////////////////////////////////////////////////////////
LittleDiracOperator LittleDiracOpPV(geom, FGrid, Coarse5d);
LittleDiracOpPV.CoarsenOperator(PVdagM, AggregatesPD);
@@ -367,7 +367,7 @@ void runMG(
//////////////////////////////////////////////////////////////////////
// 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.
// Used as the aggregation basis for Level 1->2 coarsening.
//////////////////////////////////////////////////////////////////////
std::vector<CoarseVector> psi_coarse(nbasis, Coarse5d);
for (int k = 0; k < nbasis; k++)
@@ -397,22 +397,22 @@ void runMG(
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 << "Coarse null matrix ||C - C^dag||/||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.
// 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).
// geometry 3x6x4x4 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<vTComplex>, so CComplex
// for the L1→L2 level must be iScalar<vTComplex>, not vTComplex.
// for the L1->L2 level must be iScalar<vTComplex>, not vTComplex.
typedef typename CoarseVector::vector_object CoarseSiteObj;
typedef iScalar<vTComplex> vTTComplex;
typedef GeneralCoarsenedMatrix<CoarseSiteObj,vTTComplex,NB> LittleDiracOperatorL2;
typedef DeprecatedGeneralCoarsenedMatrix<CoarseSiteObj,vTTComplex,NB> LittleDiracOperatorL2;
typedef typename LittleDiracOperatorL2::CoarseVector CoarseCoarseVector;
typedef Aggregation<CoarseSiteObj,vTTComplex,NB> SubspaceL2;
typedef MGPreconditioner<CoarseSiteObj,vTTComplex,NB> L1to2MG;
@@ -430,12 +430,12 @@ void runMG(
TrivialPrecon<CoarseCoarseVector> simpleCC;
//////////////////////////////////////////////////////////////////////
// Lüscher deflation guesser for L3PGCR.
// Luscher deflation guesser for L3PGCR.
// Step 1: project psi_coarse[k] (promoted fine null vectors) to
// CoarseCoarseVector space — these cover the zero-momentum
// 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
// 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} = <psi_cc[s]|LinOpCC|psi_cc[t]> over the
// full augmented basis and invert directly via Eigen LU.
@@ -477,13 +477,13 @@ void runMG(
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;
std::cout << GridLogMessage << "Coarse-coarse deflation matrix ||Ccc-Ccc^dag||/||Ccc|| = " << normCccmCdag/normCcc << std::endl;
}
Eigen::MatrixXcd Ccc_inv = Ccc.inverse();
LuscherGuesser<CoarseCoarseVector> CCDeflGuesser(psi_cc, Ccc_inv);
//////////////////////////////////////////////////////////////////////
// Level 2→3: coarsen LinOpCC using the RAW promoted psi_cc as aggregation
// 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
@@ -492,11 +492,11 @@ void runMG(
// so assign COPIES of psi_cc and keep psi_cc itself raw.
//
// Tensor depth deepens once more: innerProduct(CoarseCoarseSiteObj,...) returns
// iScalar<vTTComplex>, so CComplex for the L2→L3 level is iScalar<iScalar<vTComplex>>.
// iScalar<vTTComplex>, so CComplex for the L2->L3 level is iScalar<iScalar<vTComplex>>.
//////////////////////////////////////////////////////////////////////
typedef typename CoarseCoarseVector::vector_object CoarseCoarseSiteObj;
typedef iScalar<vTTComplex> vTTTComplex;
typedef GeneralCoarsenedMatrix<CoarseCoarseSiteObj,vTTTComplex,NB> LittleDiracOperatorL3;
typedef DeprecatedGeneralCoarsenedMatrix<CoarseCoarseSiteObj,vTTTComplex,NB> LittleDiracOperatorL3;
typedef typename LittleDiracOperatorL3::CoarseVector CoarseCoarseCoarseVector;
typedef Aggregation<CoarseCoarseSiteObj,vTTTComplex,NB> SubspaceL3;
typedef MGPreconditioner<CoarseCoarseSiteObj,vTTTComplex,NB> L2to3MG;
@@ -531,7 +531,7 @@ void runMG(
L4PGCR.Name("CCCouter");
//////////////////////////////////////////////////////////////////////
// Level 2→3 V-cycle: depth-2 SHIFTED smoother on LinOpCC + L4 bottom.
// 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).
//////////////////////////////////////////////////////////////////////
@@ -556,7 +556,7 @@ void runMG(
simpleCCC); // trivial guesser at the bottom
//////////////////////////////////////////////////////////////////////
// Level 3 (coarse-coarse) solve: GCR preconditioned by the L2→L3 V-cycle.
// 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.
//////////////////////////////////////////////////////////////////////
@@ -565,7 +565,7 @@ void runMG(
L3MGsolver.Name("CCouter");
//////////////////////////////////////////////////////////////////////
// Coarse-level GCR smoother for Level 1→2 V-cycle.
// 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.
//////////////////////////////////////////////////////////////////////
@@ -581,15 +581,15 @@ void runMG(
CoarseSmootherGCR.Name("Csmoother");
//////////////////////////////////////////////////////////////////////
// Level 1→2 V-cycle preconditioner.
// 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†
L3MGsolver, // coarse-coarse solve is now the L2->L3 V-cycle
CCDeflGuesser); // Luscher guesser: psi_cc C^{-1} psi_cc^dag
//////////////////////////////////////////////////////////////////////
// Standalone Level 1 two-level solve test.
@@ -624,7 +624,7 @@ void runMG(
f_src = one;
// Pre-smoother: none (TrivialPrecon); post-smoother: shifted PGCR.
// Coarse solver: L2MGsolver (PGCR preconditioned by Level 1→2 V-cycle).
// Coarse solver: L2MGsolver (PGCR preconditioned by Level 1->2 V-cycle).
TwoLevelMG ThreeLevelPrecon(AggregatesPD,
PVdagM,
simple_fine,
@@ -668,7 +668,7 @@ int main (int argc, char ** argv)
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)
// Level 1 coarse grid: block 2^4 from fine (48x48x48x96 -> 24x24x24x48, 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;
@@ -676,13 +676,13 @@ int main (int argc, char ** argv)
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).
// Level 2 coarse-coarse grid: block 2,2,3,3 from Level 1 (24x24x24x48 -> 12x12x8x16, 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.
// 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). ✓
// 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;
@@ -728,7 +728,7 @@ int main (int argc, char ** argv)
typedef PVdagMLinearOperator<MobiusFermionD,LatticeFermionD> PVdagM_t;
typedef ShiftedPVdagMLinearOperator<MobiusFermionD,LatticeFermionD> ShiftedPVdagM_t;
typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef DeprecatedGeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef LittleDiracOperator::CoarseVector CoarseVector;
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
typedef MGPreconditioner<vSpinColourVector,vTComplex,nbasis> TwoLevelMG;
examples/Example_pvdagm_5level.cc
+38 -38
View File
@@ -190,8 +190,8 @@ public:
}
};
// Lüscher deflated guesser (arXiv:0706.2298 Sec A.3) for a non-Hermitian solve.
// C_{st} = <psi[s] | LinOp | psi[t]>; guess = sum_s c_s psi[s] where c = C^{-1} psi† src.
// Luscher deflated guesser (arXiv:0706.2298 Sec A.3) for a non-Hermitian solve.
// C_{st} = <psi[s] | LinOp | psi[t]>; guess = sum_s c_s psi[s] where c = C^{-1} psi^dag src.
template<class Field>
class LuscherGuesser : public LinearFunction<Field> {
const std::vector<Field> &psi;
@@ -344,7 +344,7 @@ void runMG(
TrivialPrecon<LatticeFermionD> simple_fine;
//////////////////////////////////////////////////////////////////////
// Level 0→1: coarsen PVdagM, build LinOpCoarse
// Level 0->1: coarsen PVdagM, build LinOpCoarse
//////////////////////////////////////////////////////////////////////
LittleDiracOperator LittleDiracOpPV(geom, FGrid, Coarse5d);
LittleDiracOpPV.CoarsenOperator(PVdagM, AggregatesPD);
@@ -368,7 +368,7 @@ void runMG(
//////////////////////////////////////////////////////////////////////
// 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.
// Used as the aggregation basis for Level 1->2 coarsening.
//////////////////////////////////////////////////////////////////////
std::vector<CoarseVector> psi_coarse(nbasis, Coarse5d);
for (int k = 0; k < nbasis; k++)
@@ -398,22 +398,22 @@ void runMG(
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 << "Coarse null matrix ||C - C^dag||/||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.
// 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).
// geometry 3x6x4x4 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<vTComplex>, so CComplex
// for the L1→L2 level must be iScalar<vTComplex>, not vTComplex.
// for the L1->L2 level must be iScalar<vTComplex>, not vTComplex.
typedef typename CoarseVector::vector_object CoarseSiteObj;
typedef iScalar<vTComplex> vTTComplex;
typedef GeneralCoarsenedMatrix<CoarseSiteObj,vTTComplex,NB> LittleDiracOperatorL2;
typedef DeprecatedGeneralCoarsenedMatrix<CoarseSiteObj,vTTComplex,NB> LittleDiracOperatorL2;
typedef typename LittleDiracOperatorL2::CoarseVector CoarseCoarseVector;
typedef Aggregation<CoarseSiteObj,vTTComplex,NB> SubspaceL2;
typedef MGPreconditioner<CoarseSiteObj,vTTComplex,NB> L1to2MG;
@@ -431,12 +431,12 @@ void runMG(
TrivialPrecon<CoarseCoarseVector> simpleCC;
//////////////////////////////////////////////////////////////////////
// Lüscher deflation guesser for L3PGCR.
// Luscher deflation guesser for L3PGCR.
// Step 1: project psi_coarse[k] (promoted fine null vectors) to
// CoarseCoarseVector space — these cover the zero-momentum
// 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
// 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} = <psi_cc[s]|LinOpCC|psi_cc[t]> over the
// full augmented basis and invert directly via Eigen LU.
@@ -478,13 +478,13 @@ void runMG(
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;
std::cout << GridLogMessage << "Coarse-coarse deflation matrix ||Ccc-Ccc^dag||/||Ccc|| = " << normCccmCdag/normCcc << std::endl;
}
Eigen::MatrixXcd Ccc_inv = Ccc.inverse();
LuscherGuesser<CoarseCoarseVector> CCDeflGuesser(psi_cc, Ccc_inv);
//////////////////////////////////////////////////////////////////////
// Level 2→3: coarsen LinOpCC using the RAW promoted psi_cc as aggregation
// 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
@@ -493,11 +493,11 @@ void runMG(
// so assign COPIES of psi_cc and keep psi_cc itself raw.
//
// Tensor depth deepens once more: innerProduct(CoarseCoarseSiteObj,...) returns
// iScalar<vTTComplex>, so CComplex for the L2→L3 level is iScalar<iScalar<vTComplex>>.
// iScalar<vTTComplex>, so CComplex for the L2->L3 level is iScalar<iScalar<vTComplex>>.
//////////////////////////////////////////////////////////////////////
typedef typename CoarseCoarseVector::vector_object CoarseCoarseSiteObj;
typedef iScalar<vTTComplex> vTTTComplex;
typedef GeneralCoarsenedMatrix<CoarseCoarseSiteObj,vTTTComplex,NB> LittleDiracOperatorL3;
typedef DeprecatedGeneralCoarsenedMatrix<CoarseCoarseSiteObj,vTTTComplex,NB> LittleDiracOperatorL3;
typedef typename LittleDiracOperatorL3::CoarseVector CoarseCoarseCoarseVector;
typedef Aggregation<CoarseCoarseSiteObj,vTTTComplex,NB> SubspaceL3;
typedef MGPreconditioner<CoarseCoarseSiteObj,vTTTComplex,NB> L2to3MG;
@@ -514,7 +514,7 @@ void runMG(
TrivialPrecon<CoarseCoarseCoarseVector> simpleCCC;
//////////////////////////////////////////////////////////////////////
// Level 3→4: coarsen LinOpCCC to build the Level 5 operator, using a
// 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
@@ -526,7 +526,7 @@ void runMG(
// 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
// Tensor depth deepens once more: CComplex for the L3->L4 level is
// iScalar<vTTTComplex>. NB5 (the coarse dimension) is independent of the
// depth -- it just makes the coarsest site vector NB5-dimensional.
//////////////////////////////////////////////////////////////////////
@@ -542,10 +542,10 @@ void runMG(
// 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:
// singular value of A restricted to the span -- eig of Q^dag A^dag AQ -- NOT the
// numerical range Q^dag 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
// whitened A^dag 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.
@@ -603,7 +603,7 @@ void runMG(
typedef typename CoarseCoarseCoarseVector::vector_object CoarseCoarseCoarseSiteObj;
typedef iScalar<vTTTComplex> vTTTTComplex;
typedef GeneralCoarsenedMatrix<CoarseCoarseCoarseSiteObj,vTTTTComplex,NB5> LittleDiracOperatorL4;
typedef DeprecatedGeneralCoarsenedMatrix<CoarseCoarseCoarseSiteObj,vTTTTComplex,NB5> LittleDiracOperatorL4;
typedef typename LittleDiracOperatorL4::CoarseVector CoarseCoarseCoarseCoarseVector;
typedef Aggregation<CoarseCoarseCoarseSiteObj,vTTTTComplex,NB5> SubspaceL4;
typedef MGPreconditioner<CoarseCoarseCoarseSiteObj,vTTTTComplex,NB5> L3to4MG;
@@ -635,7 +635,7 @@ void runMG(
L5PGCR.Name("CCCCouter");
//////////////////////////////////////////////////////////////////////
// Level 3→4 V-cycle: depth-2 SHIFTED smoother on LinOpCCC + Level 5 bottom.
// 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.
//////////////////////////////////////////////////////////////////////
@@ -660,14 +660,14 @@ void runMG(
simpleCCCC); // trivial guesser at the bottom
//////////////////////////////////////////////////////////////////////
// Level 4 (coarse-coarse-coarse) solve: GCR preconditioned by the L3→L4 V-cycle.
// Level 4 (coarse-coarse-coarse) solve: GCR preconditioned by the L3->L4 V-cycle.
//////////////////////////////////////////////////////////////////////
PrecGeneralisedConjugateResidualNonHermitian<CoarseCoarseCoarseVector> 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.
// 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).
//////////////////////////////////////////////////////////////////////
@@ -688,11 +688,11 @@ void runMG(
simpleCC, // no pre-smoother
CoarseCoarseSmootherGCR, // post-smoother: depth-2 shifted GCR
LinOpCCC,
L4MGsolver, // coarse solve is now the L3→L4 V-cycle
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.
// 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.
//////////////////////////////////////////////////////////////////////
@@ -701,7 +701,7 @@ void runMG(
L3MGsolver.Name("CCouter");
//////////////////////////////////////////////////////////////////////
// Coarse-level GCR smoother for Level 1→2 V-cycle.
// 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.
//////////////////////////////////////////////////////////////////////
@@ -717,15 +717,15 @@ void runMG(
CoarseSmootherGCR.Name("Csmoother");
//////////////////////////////////////////////////////////////////////
// Level 1→2 V-cycle preconditioner.
// 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†
L3MGsolver, // coarse-coarse solve is now the L2->L3 V-cycle
CCDeflGuesser); // Luscher guesser: psi_cc C^{-1} psi_cc^dag
//////////////////////////////////////////////////////////////////////
// Standalone Level 1 two-level solve test.
@@ -760,7 +760,7 @@ void runMG(
f_src = one;
// Pre-smoother: none (TrivialPrecon); post-smoother: shifted PGCR.
// Coarse solver: L2MGsolver (PGCR preconditioned by Level 1→2 V-cycle).
// Coarse solver: L2MGsolver (PGCR preconditioned by Level 1->2 V-cycle).
TwoLevelMG ThreeLevelPrecon(AggregatesPD,
PVdagM,
simple_fine,
@@ -804,7 +804,7 @@ int main (int argc, char ** argv)
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)
// Level 1 coarse grid: block 2^4 from fine (48x48x48x96 -> 24x24x24x48, 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;
@@ -812,13 +812,13 @@ int main (int argc, char ** argv)
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).
// Level 2 coarse-coarse grid: block 2,2,3,3 from Level 1 (24x24x24x48 -> 12x12x8x16, 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.
// 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). ✓
// 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;
@@ -880,7 +880,7 @@ int main (int argc, char ** argv)
typedef PVdagMLinearOperator<MobiusFermionD,LatticeFermionD> PVdagM_t;
typedef ShiftedPVdagMLinearOperator<MobiusFermionD,LatticeFermionD> ShiftedPVdagM_t;
typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef DeprecatedGeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef LittleDiracOperator::CoarseVector CoarseVector;
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
typedef MGPreconditioner<vSpinColourVector,vTComplex,nbasis> TwoLevelMG;
examples/Example_pvdagm_census.cc
+2 -2
View File
@@ -49,7 +49,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
// sigma_min << min|lambda| : non-normal near origin
// lambda_min(H) < 0 : half-plane condition violated
//
// Requires the dagger code path in GeneralCoarsenedMatrix:
// Requires the dagger code path in DeprecatedGeneralCoarsenedMatrix:
// _Adag allocated, PopulateAdag active, _Adag exchanged, hermitian=0.
//
// Env vars:
@@ -308,7 +308,7 @@ int main (int argc, char ** argv)
typedef PVdagMLinearOperator<MobiusFermionD,LatticeFermionD> PVdagM_t;
typedef ShiftedPVdagMLinearOperator<MobiusFermionD,LatticeFermionD> ShiftedPVdagM_t;
typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef DeprecatedGeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef LittleDiracOperator::CoarseVector CoarseVector;
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
examples/Example_pvdagm_coherence.cc
+1 -1
View File
@@ -275,7 +275,7 @@ int main (int argc, char ** argv)
MobiusFermionD Dpv (Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,1.0, M5,b,c);
typedef PVdagMLinearOperator<MobiusFermionD,LatticeFermionD> PVdagM_t;
typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef DeprecatedGeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef LittleDiracOperator::CoarseVector CoarseVector;
// Grid index contraction is positional (colour/spin/lorentz order is meaningful),
// so each MG projection adds one index to the tensor nest rather than reusing a slot:
examples/Example_pvdagm_mrhs.cc
+3 -3
View File
@@ -38,7 +38,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
// precon: V-cycle -- per-RHS fine post-smoother (16-step shifted GCR),
// batched restriction (MultiRHSBlockProject / GEMM),
// ONE coarse PGCR on the 6D mrhs coarse operator
// (MultiGeneralCoarsenedMatrix, GEMM mults -- the ~10x win),
// (DeprecatedMultiGeneralCoarsenedMatrix, GEMM mults -- the ~10x win),
// batched prolongation.
//
// The coarse operator is coarsened once with the standard single-RHS
@@ -260,8 +260,8 @@ int main (int argc, char ** argv)
typedef PVdagMLinearOperator<MobiusFermionD,LatticeFermionD> PVdagM_t;
typedef ShiftedPVdagMLinearOperator<MobiusFermionD,LatticeFermionD> ShiftedPVdagM_t;
typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef MultiGeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> MrhsLittleDiracOperator;
typedef DeprecatedGeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef DeprecatedMultiGeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> MrhsLittleDiracOperator;
typedef LittleDiracOperator::CoarseVector CoarseVector;
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
examples/Example_pvdagm_mrhs_3level.cc
+5 -5
View File
@@ -23,7 +23,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
// 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.
// through GEMM (DeprecatedMultiGeneralCoarsenedMatrix) -- 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<LatticeFermionD>, MrhsPGCRNonHermitian on PVdagM,
@@ -227,16 +227,16 @@ int main (int argc, char ** argv)
typedef ShiftedPVdagMLinearOperator<MobiusFermionD,LatticeFermionD> ShiftedPVdagM_t;
// Level 1 tensor types
typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef MultiGeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> MrhsLittleDiracOperator;
typedef DeprecatedGeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef DeprecatedMultiGeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> MrhsLittleDiracOperator;
typedef LittleDiracOperator::CoarseVector CoarseVector;
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
// Level 2 tensor types (coarsening deepens the nest by one iScalar -- see CLAUDE.md)
typedef CoarseVector::vector_object CoarseSiteObj;
typedef iScalar<vTComplex> vTTComplex;
typedef GeneralCoarsenedMatrix<CoarseSiteObj,vTTComplex,nbasis> LittleDiracOperatorL2;
typedef MultiGeneralCoarsenedMatrix<CoarseSiteObj,vTTComplex,nbasis> MrhsLittleDiracOperatorL2;
typedef DeprecatedGeneralCoarsenedMatrix<CoarseSiteObj,vTTComplex,nbasis> LittleDiracOperatorL2;
typedef DeprecatedMultiGeneralCoarsenedMatrix<CoarseSiteObj,vTTComplex,nbasis> MrhsLittleDiracOperatorL2;
typedef LittleDiracOperatorL2::CoarseVector CoarseCoarseVector;
typedef Aggregation<CoarseSiteObj,vTTComplex,nbasis> SubspaceL2;
examples/Example_pvdagm_mrhs_3level_DenseCoarseMatrix.cc
+4 -4
View File
@@ -233,16 +233,16 @@ int main (int argc, char ** argv)
typedef ShiftedPVdagMLinearOperator<MobiusFermionD,LatticeFermionD> ShiftedPVdagM_t;
// Level 1 tensor types
typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef MultiGeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> MrhsLittleDiracOperator;
typedef DeprecatedGeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef DeprecatedMultiGeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> MrhsLittleDiracOperator;
typedef LittleDiracOperator::CoarseVector CoarseVector;
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
// Level 2 tensor types (coarsening deepens the nest by one iScalar)
typedef CoarseVector::vector_object CoarseSiteObj;
typedef iScalar<vTComplex> vTTComplex;
typedef GeneralCoarsenedMatrix<CoarseSiteObj,vTTComplex,nbasis> LittleDiracOperatorL2;
typedef MultiGeneralCoarsenedMatrix<CoarseSiteObj,vTTComplex,nbasis> MrhsLittleDiracOperatorL2;
typedef DeprecatedGeneralCoarsenedMatrix<CoarseSiteObj,vTTComplex,nbasis> LittleDiracOperatorL2;
typedef DeprecatedMultiGeneralCoarsenedMatrix<CoarseSiteObj,vTTComplex,nbasis> MrhsLittleDiracOperatorL2;
typedef LittleDiracOperatorL2::CoarseVector CoarseCoarseVector;
typedef Aggregation<CoarseSiteObj,vTTComplex,nbasis> SubspaceL2;
examples/Example_pvdagm_mrhs_3level_dense.cc
+5 -5
View File
@@ -237,7 +237,7 @@ public:
#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) —
// 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.
@@ -891,16 +891,16 @@ int main (int argc, char ** argv)
typedef ShiftedPVdagMLinearOperator<MobiusFermionD,LatticeFermionD> ShiftedPVdagM_t;
// Level 1 tensor types
typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef MultiGeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> MrhsLittleDiracOperator;
typedef DeprecatedGeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef DeprecatedMultiGeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> MrhsLittleDiracOperator;
typedef LittleDiracOperator::CoarseVector CoarseVector;
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
// Level 2 tensor types (coarsening deepens the nest by one iScalar)
typedef CoarseVector::vector_object CoarseSiteObj;
typedef iScalar<vTComplex> vTTComplex;
typedef GeneralCoarsenedMatrix<CoarseSiteObj,vTTComplex,nbasis> LittleDiracOperatorL2;
typedef MultiGeneralCoarsenedMatrix<CoarseSiteObj,vTTComplex,nbasis> MrhsLittleDiracOperatorL2;
typedef DeprecatedGeneralCoarsenedMatrix<CoarseSiteObj,vTTComplex,nbasis> LittleDiracOperatorL2;
typedef DeprecatedMultiGeneralCoarsenedMatrix<CoarseSiteObj,vTTComplex,nbasis> MrhsLittleDiracOperatorL2;
typedef LittleDiracOperatorL2::CoarseVector CoarseCoarseVector;
typedef Aggregation<CoarseSiteObj,vTTComplex,nbasis> SubspaceL2;
examples/Example_pvdagm_v2_3level_DenseCoarseMatrix.cc
-984
View File
@@ -1,984 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./examples/Example_pvdagm_v2_3level_DenseCoarseMatrix.cc
Copyright (C) 2026
Author: Peter Boyle <pboyle@bnl.gov>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
//
// PVdagM three level multigrid on the V2 coarse operator.
//
// STAGES ONE AND TWO: grids, types, subspace, and the L1 and L2 coarsenings.
// The dense bottom and the solves are not here yet.
//
// Differences from Example_pvdagm_mrhs_3level_DenseCoarseMatrix.cc:
//
// * The coarse space is UNVECTORISED (sComplexD). The fine space stays
// vectorised. MultiRHSBlockProject carries the mixed layout.
//
// * One operator, not two. V1 needed GeneralCoarsenedMatrix to coarsen and
// MultiGeneralCoarsenedMatrix to apply, bridged by CopyMatrix. V2 does
// both, and single versus multiRHS is SetGrid on the same object with the
// matrix elements built once.
//
// * Nrhs is unconstrained. V1 required nrhs % vComplex::Nsimd() == 0 because
// its multiRHS grid carried the SIMD in the rhs direction.
//
// * CoarsenOperator takes the subspace vectors, not an Aggregation. It block
// orthonormalises them IN PLACE -- the vectors are far too large to copy
// defensively -- so rawNull is taken first and the RAW vectors are what
// define the L2 null space. Do not insert an Orthogonalise() anywhere:
// projecting a block-orthonormal vector onto its own block-orthonormalised
// aggregation gives e_k, and the near null content is silently gone. The
// ||<psi|psi> - I||_F guard below is what catches that.
//
// Env: LATT LS MASS NBASIS(compile time) NRHS BLOCK BLOCK2 COARSEN_BATCH
// HOT_START CONFIG SUBSPACE_FILE V1_CHECK MRHS_COARSEN
//
#include <Grid/Grid.h>
#include <algorithm> // std::sort, for the runtime-environment dump in ParseEnvironment
#include <Grid/lattice/PaddedCell.h>
#include <Grid/stencil/GeneralLocalStencil.h>
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h>
#include <Grid/algorithms/multigrid/MrhsMultiGrid.h>
#include <Grid/algorithms/multigrid/DenseCoarseMatrix.h>
#include <memory>
using namespace std;
using namespace Grid;
// Compile time so it can be cut down for laptop runs: -DNBASIS=8
#ifndef NBASIS
#define NBASIS 60
#endif
RealD mass = 0.00078;
int Nrhs = 12;
int Ls = 24;
int CoarsenBatch = 9;
std::vector<int> lat_size({48,48,48,96});
// Solver tuning. PRINCIPLE (PB, 2026-08-24): the defaults ARE the current
// optimum, so an unset environment reproduces the best banked result; they
// are updated as and when a better point is found, and every change is
// dated here. Environment variables of the same names override for sweeps.
//
// Current optimum: 2026-08-24, slurm-5335492 F4, 48^3x96 Ls=24 on 288 GCDs,
// 17.2 s/RHS at Nrhs=4, 32.2 s at Nrhs=1 (exact-halo FINAL ~1e-8 pending
// the exact-outer rerun). Smoother mmax == order (full GCR history);
// PB's mmax=1 trial gave 72 vs ~60 outer iterations and was slower.
RealD FineSmootherShift = 0.1;
int FineSmootherOrder = 6;
int FineSmootherMmax = 6;
RealD CoarseSmootherShift = 0.1;
int PowerIterations = 0; // >0: power-iterate the smoother operators before the solves (spectral edge)
// Smoother implementation per level (Smoothers.h):
// gcr : the adaptive PGCR (default, as always)
// replay : run the PGCR with a coefficient recorder for the first
// PolyRecordIters outer steps, then switch to GCRReplaySmoother
// (same polynomial, no inner products) -- 1402.2585 p.13 revisited
// cheb : ChebyshevNonHermitianSmoother, 1/x on [ChebLo,ChebHi], order
// = the GCR step count of that level
std::string FineSmootherMode = "gcr";
std::string CoarseSmootherMode = "gcr";
int PolyRecordIters = 4;
int PolyRecordStart = 0; // outer step at which recording begins (0: from the first step)
std::string PolyRecordSelect = "last"; // which recorded call to replay: last|first|mean (mean is the bad one)
int PolyRefresh = 0; // >0: every PolyRefresh outer steps, one adaptive step re-records the polynomial (HDCG: every 10)
int PolyVerbose = 0; // 1: fixed-polynomial smoothers print |r_m|/|r_0| per call
RealD FineChebLo = 3.0, FineChebHi = 137.0; // from the harvested polynomial and PowerIterations edge
RealD CoarseChebLo = 8.0, CoarseChebHi = 45.0;
int CoarseSmootherNstep = 2;
int CoarseSmootherMmax = 2;
RealD CoarseSolverTol = 0.05;
int CoarseSolverOrder = 200;
int CoarseSolverMmax = 16;
RealD OuterTol = 1.0e-8;
int OuterMmax = 4;
int OuterNstep = 8;
// "It's legal to get the same answer faster, not to get a less correct
// answer." (PB, 2026-08-24)
//
// Halo-precision POLICY: reduced-precision (fp32 wire)
// halos belong in the PRECONDITIONER -- the smoother, the V-cycle's own
// residuals, and the coarsening -- and NEVER in the outer Krylov. The
// outer operator's applications define what "converged" means; making them
// sloppy turns the stopping criterion into a statement about the wrong
// operator (measured: solver stops at computed 9.8e-9 while the true
// residual is 3.4e-8). Exactness costs one exact fine matvec per outer
// iteration, a few percent of the solve.
//
// Stencil::SloppyComms is a free runtime setter (Stencil.h:303), so the
// policy is implemented by SCOPED toggling: SetFineSloppy(1) on entering
// the preconditioner / coarsening, SetFineSloppy(0) on leaving. The
// operators default to EXACT. FineSloppyComms therefore now means
// "sloppy inside the preconditioner"; =0 makes everything exact.
int FineSloppyComms = 1;
// SmootherCoeffLog=1 : print the GCR step lengths a_k and orthogonalisation
// coefficients b_kj of BOTH smoothers every call -- the harvest for a fixed
// polynomial smoother (stable coefficients => stationary p(A), no reductions).
int SmootherCoeffLog = 0;
std::function<void(int)> SetFineSloppy = [](int){};
void ParseEnvironment(void)
{
if(getenv("MASS")) mass = atof(getenv("MASS"));
if(getenv("NRHS")) Nrhs = atoi(getenv("NRHS"));
if(getenv("LS")) Ls = atoi(getenv("LS"));
if(getenv("COARSEN_BATCH")) CoarsenBatch= atoi(getenv("COARSEN_BATCH"));
if(getenv("FineSmootherShift")) FineSmootherShift = atof(getenv("FineSmootherShift"));
if(getenv("FineSmootherOrder")) FineSmootherOrder = atoi(getenv("FineSmootherOrder"));
if(getenv("FineSmootherMmax")) FineSmootherMmax = atoi(getenv("FineSmootherMmax"));
if(getenv("CoarseSmootherShift"))CoarseSmootherShift= atof(getenv("CoarseSmootherShift"));
if(getenv("PowerIterations")) PowerIterations = atoi(getenv("PowerIterations"));
if(getenv("FineSmootherMode")) FineSmootherMode = getenv("FineSmootherMode");
if(getenv("CoarseSmootherMode")) CoarseSmootherMode = getenv("CoarseSmootherMode");
if(getenv("PolyRecordIters")) PolyRecordIters = atoi(getenv("PolyRecordIters"));
if(getenv("PolyRecordStart")) PolyRecordStart = atoi(getenv("PolyRecordStart"));
if(getenv("PolyRecordSelect")) PolyRecordSelect = getenv("PolyRecordSelect");
if(getenv("PolyRefresh")) PolyRefresh = atoi(getenv("PolyRefresh"));
if(getenv("PolyVerbose")) PolyVerbose = atoi(getenv("PolyVerbose"));
if(getenv("FineChebLo")) FineChebLo = atof(getenv("FineChebLo"));
if(getenv("FineChebHi")) FineChebHi = atof(getenv("FineChebHi"));
if(getenv("CoarseChebLo")) CoarseChebLo = atof(getenv("CoarseChebLo"));
if(getenv("CoarseChebHi")) CoarseChebHi = atof(getenv("CoarseChebHi"));
if(getenv("CoarseSmootherNstep"))CoarseSmootherNstep= atoi(getenv("CoarseSmootherNstep"));
if(getenv("CoarseSmootherMmax")) CoarseSmootherMmax = atoi(getenv("CoarseSmootherMmax"));
if(getenv("CoarseSolverTol")) CoarseSolverTol = atof(getenv("CoarseSolverTol"));
if(getenv("CoarseSolverOrder")) CoarseSolverOrder = atoi(getenv("CoarseSolverOrder"));
if(getenv("CoarseSolverMmax")) CoarseSolverMmax = atoi(getenv("CoarseSolverMmax"));
if(getenv("OuterTol")) OuterTol = atof(getenv("OuterTol"));
if(getenv("OuterMmax")) OuterMmax = atoi(getenv("OuterMmax"));
if(getenv("FineSloppyComms")) FineSloppyComms = atoi(getenv("FineSloppyComms"));
if(getenv("SmootherCoeffLog")) SmootherCoeffLog = atoi(getenv("SmootherCoeffLog"));
if(getenv("OuterNstep")) OuterNstep = atoi(getenv("OuterNstep"));
if(getenv("LATT")){
Coordinate l;
GridCmdOptionIntVector(std::string(getenv("LATT")),l);
GRID_ASSERT(l.size()==4);
for(int d=0;d<4;d++) lat_size[d]=l[d];
}
std::cout << GridLogMessage << "PARAM: LATT "
<< lat_size[0]<<"."<<lat_size[1]<<"."<<lat_size[2]<<"."<<lat_size[3] << std::endl;
std::cout << GridLogMessage << "PARAM: LS " << Ls << std::endl;
std::cout << GridLogMessage << "PARAM: MASS " << mass << std::endl;
std::cout << GridLogMessage << "PARAM: NBASIS " << NBASIS << std::endl;
std::cout << GridLogMessage << "PARAM: NRHS " << Nrhs << std::endl;
std::cout << GridLogMessage << "PARAM: COARSEN_BATCH " << CoarsenBatch << std::endl;
// EVERY knob this programme parsed, so a log identifies its own run (2026-08-28:
// four jobs in the queue differing in FineSloppyComms/NRHS/mmax and none of it
// printed). Interim until the Serializable parameter struct replaces all of this.
auto P = [](const char *n, auto v){ std::cout << GridLogMessage << "PARAM: " << std::left << std::setw(20) << n << std::right << " " << v << std::endl; };
P("FineSmootherShift",FineSmootherShift); P("FineSmootherOrder",FineSmootherOrder); P("FineSmootherMmax",FineSmootherMmax);
P("FineSmootherMode",FineSmootherMode); P("FineChebLo",FineChebLo); P("FineChebHi",FineChebHi);
P("FineSloppyComms",FineSloppyComms);
P("CoarseSmootherShift",CoarseSmootherShift); P("CoarseSmootherNstep",CoarseSmootherNstep); P("CoarseSmootherMmax",CoarseSmootherMmax);
P("CoarseSmootherMode",CoarseSmootherMode); P("CoarseChebLo",CoarseChebLo); P("CoarseChebHi",CoarseChebHi);
P("CoarseSolverTol",CoarseSolverTol); P("CoarseSolverOrder",CoarseSolverOrder); P("CoarseSolverMmax",CoarseSolverMmax);
P("OuterTol",OuterTol); P("OuterMmax",OuterMmax); P("OuterNstep",OuterNstep);
P("PowerIterations",PowerIterations); P("PolyRecordIters",PolyRecordIters); P("PolyRecordStart",PolyRecordStart);
P("PolyRecordSelect",PolyRecordSelect); P("PolyRefresh",PolyRefresh); P("PolyVerbose",PolyVerbose);
P("SmootherCoeffLog",SmootherCoeffLog);
// library-side knobs, read straight from the environment as the library will
const char *envs[] = {"BLOCK","BLOCK2","SUBSPACE_FILE","CONFIG","HOT_START","MRHS_COARSEN","V1_CHECK",
"DENSE_CC","DENSE_SCHUR","DENSE_SCHUR2D","DENSE_DEVICE_SUM","DENSE_SPLITK","DENSE_NB",
"SCHUR2D_LEAF_SPAN","SCHUR2D_LEAF_LU","SCHUR2D_PROBE","SUMMA_HANDSHAKE","DENSE_APPLY_PROFILE","SLAB_FILE",
"SOLVE_SRHS"};
for(auto e : envs) P(e, getenv(e) ? std::string(getenv(e)) : std::string("(unset)"));
////////////////////////////////////////////////////////////////////////////////////////
// The COMPLETE fabric/runtime environment, scanned from environ rather than a curated
// list. 2026-08-30: we could not tell from a failing log whether FI_CXI_ATS was set,
// because it was not in the list -- and a module or the submitting shell can set
// anything. A log must describe its own run without reference to the job script.
////////////////////////////////////////////////////////////////////////////////////////
{
extern char **environ;
const char *prefixes[] = {"FI_","OFI_","MPICH_","MPIR_","PMI_","CXI_","HSA_","HIP_","ROCR_",
"ROCM_","AMD_","GPU_","NUMA_","OMP_","GOMP_","GRID_","CRAY_","LIBFABRIC"};
std::vector<std::string> hits;
for(char **e = environ; e && *e; e++){
std::string s(*e);
for(auto p : prefixes){
if( s.compare(0,strlen(p),p)==0 ){
if(s.size()>200) s = s.substr(0,197)+"..."; // paths can be enormous
hits.push_back(s);
break;
}
}
}
std::sort(hits.begin(),hits.end());
std::cout << GridLogMessage << "PARAM: ---- runtime environment: "<<hits.size()<<" variables set ----"<<std::endl;
for(auto &s : hits) std::cout << GridLogMessage << "PARAM: ENV " << s << std::endl;
// Variables whose ABSENCE is as informative as their value (defaults bite: an unset
// FI_MR_CACHE_MONITOR means libfabric's defective memhooks, see systems/WorkArounds.txt).
const char *notable[] = {"FI_MR_CACHE_MONITOR","FI_MR_CACHE_MAX_COUNT","FI_MR_ROCR_CACHE_MONITOR_ENABLED",
"FI_CXI_ATS","FI_CXI_RDZV_THRESHOLD","FI_CXI_DISABLE_HMEM_DEV_REGISTER",
"FI_HMEM_ROCR_USE_DMABUF","MPICH_GPU_SUPPORT_ENABLED","MPICH_OFI_NIC_POLICY",
"MPICH_SMP_SINGLE_COPY_MODE","OMP_NUM_THREADS","GRID_ALLOC_NCACHE_LARGE"};
for(auto n : notable) if(!getenv(n))
std::cout << GridLogMessage << "PARAM: ENV " << n << " (UNSET - provider/runtime default applies)" << std::endl;
}
}
template <class Field>
void saveSubspace(std::vector<Field> &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 <class Field>
void loadSubspace(std::vector<Field> &subspace, std::string const fname){
#ifdef HAVE_LIME
Grid::emptyUserRecord record;
Grid::ScidacReader SR;
SR.open(fname);
for (int k = 0; k < (int)subspace.size(); k++) SR.readScidacFieldRecord(subspace[k], record);
SR.close();
#endif
}
//////////////////////////////////////////////////////////////////////
// A = PV^dag M (non-Hermitian)
//////////////////////////////////////////////////////////////////////
template<class Matrix,class Field>
class PVdagMLinearOperator : public LinearOperatorBase<Field> {
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<Field> &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); }
};
//////////////////////////////////////////////////////////////////////
// ||<v|v> - I||_F over a set of coarse vectors. Small means the raw near null
// content survived the projection; see GramGuard for where a leak lands.
//////////////////////////////////////////////////////////////////////
template<class CoarseField>
RealD GramDefect(std::vector<CoarseField> &v)
{
RealD s2=0.0;
for(int i=0;i<(int)v.size();i++){
for(int j=0;j<(int)v.size();j++){
ComplexD sij=TensorRemove(innerProduct(v[i],v[j]));
ComplexD d=sij-(i==j?ComplexD(1.0):ComplexD(0.0));
s2+=real(d)*real(d)+imag(d)*imag(d);
}
}
return std::sqrt(s2);
}
// On a leak every image collapses to the block unit e_k, the Gram becomes
// N*I, and the defect lands at (N-1)*sqrt(nbasis) -- orders above the ~0.2
// of a content preserving projection. Trip well below that so a mis-set
// threshold costs a log line rather than the run.
template<class CoarseField>
void GramGuard(const std::string &name,std::vector<CoarseField> &v,GridBase *grid)
{
RealD defect = GramDefect(v);
RealD N = (RealD)grid->gSites();
RealD leak = (N-1.0)*std::sqrt((RealD)v.size());
RealD trip = std::sqrt(N);
std::cout << GridLogMessage << "GUARD: ||<"<<name<<"|"<<name<<"> - I||_F = " << defect
<< " (e_k leak would be " << leak << ", trip at " << trip << ")" << std::endl;
GRID_ASSERT( defect < trip );
}
//////////////////////////////////////////////////////////////////////
// Shifted variants for the smoothers
//////////////////////////////////////////////////////////////////////
template<class Matrix,class Field>
class ShiftedPVdagMLinearOperator : public LinearOperatorBase<Field> {
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<Field> &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 Field>
class ShiftedLinearOperator : public LinearOperatorBase<Field> {
LinearOperatorBase<Field> &_Op; RealD shift;
public:
ShiftedLinearOperator(RealD _shift, LinearOperatorBase<Field> &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<Field> &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); }
};
//////////////////////////////////////////////////////////////////////
// Power iteration on a (non-Hermitian) operator: the spectral edge the
// smoother polynomial must not exceed. Reports
// step 0 : |A v|/|v| on a RANDOM unit v -- a one-sample lower bound on
// sigma_max(A). If this and the converged value agree, the
// operator is near-normal and the spectral picture (R_m(lambda)
// on the spectrum) is trustworthy; if not, the field of values
// sets the safe interval and the spectrum understates it.
// step k : |A v_k|/|v_k| -> |lambda_max| as v_k -> the dominant
// eigenvector; the complex Rayleigh quotient <v,Av> gives its
// phase (real => on the axis). A non-converging oscillation
// means a complex-conjugate pair of equal modulus at the top.
// Uses Op(), not HermOp(): this is the operator the smoother sees.
//////////////////////////////////////////////////////////////////////
template<class Field>
void PowerIteration(const std::string &name, LinearOperatorBase<Field> &Op, GridBase *grid, int iters)
{
GRID_TRACE("PowerIteration");
GridParallelRNG RNG(grid); RNG.SeedFixedIntegers(std::vector<int>({7,11,13,17}));
Field v(grid), Av(grid);
gaussian(RNG,v);
RealD nv = std::sqrt(norm2(v)); v = v*(1.0/nv);
RealD ratio=0.0, ratio0=0.0; ComplexD rq(0.0);
for(int i=0;i<iters;i++){
Op.Op(v,Av);
RealD nAv = std::sqrt(norm2(Av));
rq = innerProduct(v,Av); // v is unit
ratio = nAv;
if ( i==0 ) {
ratio0 = ratio;
std::cout << GridLogMessage << "PowerIteration " << name << " step 0 (random v): |Av|/|v| = " << ratio
<< " [lower bound on sigma_max]" << std::endl;
}
if ( (i%10==0) || (i==iters-1) )
std::cout << GridLogMessage << "PowerIteration " << name << " step " << i << " |Av|/|v| = " << ratio
<< " Rayleigh = (" << real(rq) << "," << imag(rq) << ")" << std::endl;
v = Av*(1.0/nAv);
}
std::cout << GridLogMessage << "PowerIteration " << name << " SUMMARY: |lambda_max| ~ " << ratio
<< " Rayleigh (" << real(rq) << "," << imag(rq) << ")"
<< " phase " << std::atan2(imag(rq),real(rq)) << " rad"
<< " step-0 ratio / converged = " << ratio0/ratio
<< (ratio0/ratio > 1.2 ? " ** non-normal: sigma_max well above |lambda_max| **" : " (near-normal)")
<< std::endl;
}
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
ParseEnvironment();
RealD M5=1.8, b=1.5, c=0.5;
const int nbasis=NBASIS;
const int nrhs=Nrhs;
const int batch=CoarsenBatch;
Coordinate mpi = GridDefaultMpi();
Coordinate fsimd= GridDefaultSimd(Nd,vComplex::Nsimd());
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(lat_size,fsimd,mpi);
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
// Level 1 blocking (default 2^4)
Coordinate clatt = lat_size;
Coordinate Block({2,2,3,3}); // L1->L2 blocking; banked optimum 2026-08-24 (env BLOCK overrides)
if ( getenv("BLOCK") ){ GridCmdOptionIntVector(std::string(getenv("BLOCK")),Block); GRID_ASSERT(Block.size()==4); }
for(int d=0;d<4;d++){ GRID_ASSERT(lat_size[d]%Block[d]==0); clatt[d]=lat_size[d]/Block[d]; }
std::cout << GridLogMessage << "Block " << Block << " coarse lattice " << clatt << std::endl;
//////////////////////////////////////////////////////////////////////
// The coarse space is unvectorised. The 5D coarse grid is built here
// rather than through SpaceTimeGrid so the SIMD layout is ours.
//////////////////////////////////////////////////////////////////////
Coordinate c5latt({1,clatt[0],clatt[1],clatt[2],clatt[3]});
Coordinate c5simd({1,1,1,1,1});
Coordinate c5mpi ({1,mpi[0],mpi[1],mpi[2],mpi[3]});
GridCartesian *Coarse5d = new GridCartesian(c5latt,c5simd,c5mpi);
// 6D coarse multiRHS grid: rhs is dim 0, undistributed and unvectorised.
// No divisibility constraint on nrhs, unlike V1.
Coordinate cmlatt({nrhs,1,clatt[0],clatt[1],clatt[2],clatt[3]});
Coordinate cmsimd({1,1,1,1,1,1});
Coordinate cmmpi ({1,1,mpi[0],mpi[1],mpi[2],mpi[3]});
GridCartesian *CoarseMrhs = new GridCartesian(cmlatt,cmsimd,cmmpi);
// 6D coarse grid at the coarsening batch, used only while CoarsenOperator
// runs. The matrix elements survive the change back to nrhs.
Coordinate cblatt({batch,1,clatt[0],clatt[1],clatt[2],clatt[3]});
GridCartesian *CoarseBatch = new GridCartesian(cblatt,cmsimd,cmmpi);
// 6D fine grid carrying the coarsening batch: fine SIMD layout preserved
Coordinate fmlatt({batch,Ls,lat_size[0],lat_size[1],lat_size[2],lat_size[3]});
Coordinate fmsimd({1,1,fsimd[0],fsimd[1],fsimd[2],fsimd[3]});
Coordinate fmmpi ({1,1,mpi[0],mpi[1],mpi[2],mpi[3]});
GridCartesian *FineMrhs = new GridCartesian(fmlatt,fmsimd,fmmpi);
std::cout << GridLogMessage << "Nsimd fine " << FGrid->Nsimd()
<< " coarse " << Coarse5d->Nsimd() << std::endl;
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers({1,2,3,4});
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers({5,6,7,8});
//////////////////////////////////////////////////////////////////////
// Gauge field
//////////////////////////////////////////////////////////////////////
LatticeGaugeField Umu(UGrid);
if ( getenv("HOT_START") ) {
std::cout << GridLogMessage << "Hot start gauge field" << std::endl;
SU<Nc>::HotConfiguration(RNG4,Umu);
} else {
std::string file("/ccs/home/poare/ckpoint_lat.1000");
if ( getenv("CONFIG") ) file = std::string(getenv("CONFIG"));
std::cout << GridLogMessage << "Reading gauge field " << file << std::endl;
FieldMetaData header;
NerscIO::readConfiguration(Umu,header,file);
}
MobiusFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c);
MobiusFermionD Dpv (Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,1.0, M5,b,c);
// PVdagM and ShiftedPVdagM are thin wrappers over these same two objects,
// so this one callback controls every fine halo in the program. Default
// EXACT; the preconditioner and the coarsening turn sloppiness on for
// their own scope only (policy note at FineSloppyComms).
SetFineSloppy = [&Ddwf,&Dpv](int sloppy){
Ddwf.SloppyComms(sloppy);
Dpv .SloppyComms(sloppy);
};
SetFineSloppy(0);
std::cout << GridLogMessage << "Fine halo policy: preconditioner+coarsening "
<< (FineSloppyComms ? "SLOPPY (fp32 wire)" : "exact")
<< ", outer Krylov EXACT" << std::endl;
typedef PVdagMLinearOperator<MobiusFermionD,LatticeFermionD> PVdagM_t;
typedef ShiftedPVdagMLinearOperator<MobiusFermionD,LatticeFermionD> ShiftedPVdagM_t;
PVdagM_t PVdagM(Ddwf,Dpv);
//////////////////////////////////////////////////////////////////////
// Level 1 types: unvectorised coarse scalar
//////////////////////////////////////////////////////////////////////
typedef sTComplexD CComplexS;
typedef MultiGeneralCoarsenedOperatorV2<vSpinColourVector,CComplexS,nbasis> CoarseOperator;
typedef CoarseOperator::CoarseVector CoarseVector;
typedef Aggregation<vSpinColourVector,CComplexS,nbasis> Subspace;
NextToNearestStencilGeometry5D geom(Coarse5d);
//////////////////////////////////////////////////////////////////////
// Subspace: load RAW (no Orthogonalise!), or generate.
//
// The Aggregation is scaffolding for CreateSubspaceGCR only. That runs
// entirely on the fine grid and ends in GlobalOrthonormalise, which is a
// whole-lattice Gram-Schmidt, so the coarse grid it holds is never
// dereferenced and may be the unvectorised one.
//////////////////////////////////////////////////////////////////////
std::string subspace_file = "subspace_nb" + std::to_string(nbasis) + ".scidac";
if ( getenv("SUBSPACE_FILE") ) subspace_file = std::string(getenv("SUBSPACE_FILE"));
uint64_t file_exists=0;
if ( UGrid->IsBoss() ){ std::ifstream f(subspace_file); file_exists=f.good()?1:0; }
UGrid->GlobalSum(file_exists);
const int cb=0;
Subspace AggregatesGCR(Coarse5d,FGrid,cb);
if ( file_exists ){
std::cout << GridLogMessage << "*** Loading subspace from disk (kept RAW) ***" << std::endl;
loadSubspace(AggregatesGCR.subspace, subspace_file);
} else {
std::cout << GridLogMessage << "*** GCR subspace generation ***" << std::endl;
AggregatesGCR.CreateSubspaceGCR(RNG5,PVdagM,nbasis);
saveSubspace(AggregatesGCR.subspace, subspace_file);
}
// RAW copy BEFORE CoarsenOperator block-orthonormalises in place.
std::vector<LatticeFermionD> rawNull(nbasis,FGrid);
for(int k=0;k<nbasis;k++) rawNull[k]=AggregatesGCR.subspace[k];
//////////////////////////////////////////////////////////////////////
// L1 coarsening. The fine operator is single RHS, so it is promoted to
// the 6D batch grid; a natively multiRHS fine operator would substitute
// here with no other change.
//////////////////////////////////////////////////////////////////////
CoarseOperator CoarseOpPV(geom,Coarse5d);
CoarseOpPV.SetGrid(CoarseBatch);
std::cout << GridLogMessage << "*** L1 CoarsenOperator, batch "<<batch<<" ***" << std::endl;
SetFineSloppy(FineSloppyComms); // coarsening builds the PRECONDITIONER
if ( getenv("MRHS_COARSEN") ) {
// Promote the single RHS operator and pack the batch: costs an
// ExtractSlice/InsertSlice pair per rhs. Here for the A/B only; this is
// the path a natively multiRHS fine operator would take.
MrhsPromotedOperator<LatticeFermionD> MrhsPVdagM(PVdagM,FGrid,batch);
CoarseOpPV.CoarsenOperator(MrhsPVdagM,FineMrhs,AggregatesGCR.subspace,Coarse5d);
} else {
// PVdagM is single RHS: apply it directly, batch on the coarse side.
CoarseOpPV.CoarsenOperator(PVdagM,AggregatesGCR.subspace,Coarse5d,batch);
}
SetFineSloppy(0);
// (moved here 2026-08-28: it needs AggregatesGCR.subspace, which is freed right after
// the projector imports it below -- the 60 fine vectors are 20 GB of host memory per
// rank and 60 LRU-eligible device fields competing with the solver's working set)
//////////////////////////////////////////////////////////////////////
// Optional cross check of the coarse matrix elements against the V1
// path, which needs a vectorised coarse space. Block Gram-Schmidt is
// idempotent, so V1 may re-orthonormalise the same vectors in place
// without a second copy of the subspace.
//////////////////////////////////////////////////////////////////////
if ( getenv("V1_CHECK") ) {
typedef GeneralCoarsenedMatrix <vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef MultiGeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> MrhsLittleDiracOperator;
typedef Aggregation <vSpinColourVector,vTComplex,nbasis> SubspaceV;
Coordinate v5latt({1,clatt[0],clatt[1],clatt[2],clatt[3]});
Coordinate v5simd({1,fsimd[0],fsimd[1],fsimd[2],fsimd[3]});
GridCartesian *Coarse5dV = new GridCartesian(v5latt,v5simd,c5mpi);
int nrhs_v1 = vComplex::Nsimd();
Coordinate vmlatt({nrhs_v1,1,clatt[0],clatt[1],clatt[2],clatt[3]});
Coordinate vmsimd({vComplex::Nsimd(),1,1,1,1,1});
GridCartesian *CoarseMrhsV = new GridCartesian(vmlatt,vmsimd,cmmpi);
NextToNearestStencilGeometry5D geomV(Coarse5dV);
SubspaceV AggV(Coarse5dV,FGrid,cb);
for(int k=0;k<nbasis;k++) AggV.subspace[k]=AggregatesGCR.subspace[k];
LittleDiracOperator LittleDiracOpPV(geomV,FGrid,Coarse5dV);
std::cout << GridLogMessage << "*** V1 CoarsenOperator (cross check) ***" << std::endl;
SetFineSloppy(FineSloppyComms);
LittleDiracOpPV.CoarsenOperator(PVdagM,AggV);
SetFineSloppy(0);
MrhsLittleDiracOperator mrhsV1(geomV,CoarseMrhsV);
mrhsV1.CopyMatrix(LittleDiracOpPV);
// BLAS_A is written by GridtoBLAS in lSite order and both sides carry
// the same scalar_object, so the two are directly comparable.
typedef MrhsLittleDiracOperator::calcMatrix calcMatrix;
int npoint = geom.npoint;
RealD num=0.0, den=0.0;
for(int p=0;p<npoint;p++){
int64_t sites = mrhsV1.BLAS_A[p].size();
GRID_ASSERT(sites == (int64_t)CoarseOpPV.BLAS_A[p].size());
std::vector<calcMatrix> h1(sites),h2(sites);
acceleratorCopyFromDevice(&mrhsV1.BLAS_A[p][0], &h1[0],sites*sizeof(calcMatrix));
acceleratorCopyFromDevice(&CoarseOpPV.BLAS_A[p][0],&h2[0],sites*sizeof(calcMatrix));
ComplexD *w1=(ComplexD *)&h1[0];
ComplexD *w2=(ComplexD *)&h2[0];
int64_t words = sites*sizeof(calcMatrix)/sizeof(ComplexD);
for(int64_t i=0;i<words;i++){
ComplexD d=w1[i]-w2[i];
num += real(d)*real(d)+imag(d)*imag(d);
den += real(w1[i])*real(w1[i])+imag(w1[i])*imag(w1[i]);
}
}
std::cout << GridLogMessage << "V1_CHECK: |A_V1|^2 = " << den << std::endl;
std::cout << GridLogMessage << "V1_CHECK: |A_V1 - A_V2|^2 / |A_V1|^2 = " << num/den << std::endl;
GRID_ASSERT( den > 0.0 );
GRID_ASSERT( num/den < 1.0e-18 );
}
// Stay on the batch grid: the L2 coarsening drives this operator at the
// batch. It is switched to the solve Nrhs once L2 is built.
//////////////////////////////////////////////////////////////////////
// psi_coarse = P^dag (RAW fine null) -> Galerkin images that carry the
// near null content, and are free: A_c (P psi) = P A psi.
//////////////////////////////////////////////////////////////////////
MultiRHSBlockProject<LatticeFermionD> MrhsProjector;
MrhsProjector.Allocate(nbasis,FGrid,Coarse5d);
MrhsProjector.ImportBasis(AggregatesGCR.subspace); // block orthonormal basis
// The basis now lives in the projector's BLAS_V (device); nothing after this
// point reads AggregatesGCR.subspace (V1_CHECK moved above). Free it: 60 fine
// fields = 20 GB host per rank, and 60 LRU-eligible device copies.
AggregatesGCR.subspace.clear(); AggregatesGCR.subspace.shrink_to_fit();
std::vector<CoarseVector> psi_coarse(nbasis,Coarse5d);
MrhsProjector.blockProject(rawNull,psi_coarse); // RAW vectors in
rawNull.clear(); rawNull.shrink_to_fit();
GramGuard("psi_coarse",psi_coarse,Coarse5d);
//////////////////////////////////////////////////////////////////////
// STAGE TWO: L2 -> L3.
//
// The fine operator here is V2 at L1, which is natively multiRHS, so the
// multiRHS driver applies with no promotion adapter: its D+1 grid IS the
// batch grid the L1 operator is currently set to.
//////////////////////////////////////////////////////////////////////
Coordinate cclatt = clatt;
Coordinate Block2({4,4,2,4}); // L2->L3 blocking; banked optimum 2026-08-24 (env BLOCK2 overrides)
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;
Coordinate cc5latt({1,cclatt[0],cclatt[1],cclatt[2],cclatt[3]});
GridCartesian *CoarseCoarse5d = new GridCartesian(cc5latt,c5simd,c5mpi);
Coordinate ccmlatt({nrhs,1,cclatt[0],cclatt[1],cclatt[2],cclatt[3]});
GridCartesian *CoarseCoarseMrhs = new GridCartesian(ccmlatt,cmsimd,cmmpi);
Coordinate ccblatt({batch,1,cclatt[0],cclatt[1],cclatt[2],cclatt[3]});
GridCartesian *CoarseCoarseBatch = new GridCartesian(ccblatt,cmsimd,cmmpi);
// Coarsening deepens the tensor nest by one iScalar
typedef CoarseVector::vector_object CoarseSiteObj;
typedef iScalar<CComplexS> CComplexS2;
typedef MultiGeneralCoarsenedOperatorV2<CoarseSiteObj,CComplexS2,nbasis> CoarseCoarseOperator;
typedef CoarseCoarseOperator::CoarseVector CoarseCoarseVector;
NextToNearestStencilGeometry5D geom2(CoarseCoarse5d);
// RAW copy of the coarse null vectors, for the same reason as rawNull:
// the L2 CoarsenOperator block-orthonormalises its subspace in place, and
// the L3 basis must be defined by the vectors that still carry content.
std::vector<CoarseVector> rawPsi(nbasis,Coarse5d);
for(int k=0;k<nbasis;k++) rawPsi[k]=psi_coarse[k];
CoarseCoarseOperator CoarseOpL2(geom2,CoarseCoarse5d);
CoarseOpL2.SetGrid(CoarseCoarseBatch);
NonHermitianLinearOperator<CoarseOperator,CoarseVector> LinOpCoarse(CoarseOpPV);
std::cout << GridLogMessage << "*** L2 CoarsenOperator, batch "<<batch<<" ***" << std::endl;
CoarseOpL2.CoarsenOperator(LinOpCoarse,CoarseBatch,psi_coarse,CoarseCoarse5d);
//////////////////////////////////////////////////////////////////////
// Both operators to the solve Nrhs. The matrix elements are Nrhs
// independent and survive the change.
//////////////////////////////////////////////////////////////////////
CoarseOpPV.SetGrid(CoarseMrhs);
CoarseOpL2.SetGrid(CoarseCoarseMrhs);
std::cout << GridLogMessage << "L1 operator at Nrhs " << CoarseOpPV.Nrhs()
<< ", L2 operator at Nrhs " << CoarseOpL2.Nrhs() << std::endl;
//////////////////////////////////////////////////////////////////////
// psi_cc from the RAW coarse null vectors, and the same guard
//////////////////////////////////////////////////////////////////////
MultiRHSBlockProject<CoarseVector> MrhsProjectorL2;
MrhsProjectorL2.Allocate(nbasis,Coarse5d,CoarseCoarse5d);
MrhsProjectorL2.ImportBasis(psi_coarse); // block orthonormal basis
{
std::vector<CoarseCoarseVector> psi_cc(nbasis,CoarseCoarse5d);
MrhsProjectorL2.blockProject(rawPsi,psi_cc); // RAW vectors in
GramGuard("psi_cc",psi_cc,CoarseCoarse5d);
}
rawPsi.clear(); rawPsi.shrink_to_fit();
//////////////////////////////////////////////////////////////////////
// Both coarse operators apply on their solve grids
//////////////////////////////////////////////////////////////////////
{
GridParallelRNG cRNG(Coarse5d); cRNG.SeedFixedIntegers({3,4,5,6});
CoarseVector cin(CoarseMrhs), cout_(CoarseMrhs);
random(cRNG,cin);
CoarseOpPV.M(cin,cout_);
std::cout << GridLogMessage << "L1 apply |in|^2 = " << norm2(cin)
<< " |M in|^2 = " << norm2(cout_) << std::endl;
GRID_ASSERT( norm2(cout_) > 0.0 );
GridParallelRNG ccRNG(CoarseCoarse5d); ccRNG.SeedFixedIntegers({7,8,9,10});
CoarseCoarseVector ccin(CoarseCoarseMrhs), ccout(CoarseCoarseMrhs);
random(ccRNG,ccin);
CoarseOpL2.M(ccin,ccout);
std::cout << GridLogMessage << "L2 apply |in|^2 = " << norm2(ccin)
<< " |M in|^2 = " << norm2(ccout) << std::endl;
GRID_ASSERT( norm2(ccout) > 0.0 );
}
//////////////////////////////////////////////////////////////////////
// STAGE THREE (part one): the dense bottom on L2.
//
// DenseCoarseMatrix is bilingual: it takes the elements through
// Geometry()/ExtractMatrix(), so the V2 operator serves directly. It does
// detect that a multiRHS op cannot apply on the D dimensional grid and
// skips its own certificate and VERIFY, so the equivalent check is done
// here instead, driving the L2 operator at Nrhs 1 through a slice.
//////////////////////////////////////////////////////////////////////
typedef DenseCoarseMatrix<CComplexS2,nbasis> DenseCC_t;
std::unique_ptr<DenseCC_t> DenseCC;
if ( getenv("DENSE_CC")==nullptr || atoi(getenv("DENSE_CC")) ) {
std::cout << GridLogMessage << "*** L3 dense bottom: import from the V2 L2 operator ***" << std::endl;
DenseCC.reset(new DenseCC_t(CoarseCoarse5d));
DenseCC->Import(CoarseOpL2);
////////////////////////////////////////////////////////////////////
// ||A Ainv x - x|| / ||x||, the check Import could not run itself
////////////////////////////////////////////////////////////////////
Coordinate cc1latt({1,1,cclatt[0],cclatt[1],cclatt[2],cclatt[3]});
GridCartesian *CoarseCoarseOne = new GridCartesian(cc1latt,cmsimd,cmmpi);
CoarseOpL2.SetGrid(CoarseCoarseOne);
CoarseCoarseVector x(CoarseCoarse5d),y(CoarseCoarse5d),z(CoarseCoarse5d);
GridParallelRNG dRNG(CoarseCoarse5d); dRNG.SeedFixedIntegers({11,12,13,14});
random(dRNG,x);
(*DenseCC)(x,y); // y = Ainv x
CoarseCoarseVector y1(CoarseCoarseOne),z1(CoarseCoarseOne);
InsertSliceFast(y,y1,0,0);
CoarseOpL2.M(y1,z1); // z = A y
ExtractSliceFast(z,z1,0,0);
z = z - x;
RealD rel = std::sqrt(norm2(z)/norm2(x));
std::cout << GridLogMessage << "L3 dense: ||A Ainv x - x||/||x|| = " << rel << std::endl;
GRID_ASSERT( rel < 1.0e-2 );
CoarseOpL2.SetGrid(CoarseCoarseMrhs);
delete CoarseCoarseOne;
}
//////////////////////////////////////////////////////////////////////
// STAGE THREE (part two): the solves.
//
// Both operators are driven from the SAME objects at whatever Nrhs is
// asked for -- the matrix elements were built once and survive SetGrid --
// so single RHS and multiRHS are the same code path with a different grid.
//////////////////////////////////////////////////////////////////////
GRID_ASSERT(DenseCC != nullptr); // the PGCR bottom is not ported yet
typedef PrecGeneralisedConjugateResidualNonHermitian<LatticeFermionD> FineSmoother_t;
ShiftedPVdagM_t ShiftedPVdagM(FineSmootherShift,Ddwf,Dpv);
TrivialPrecon<LatticeFermionD> simple_fine;
TrivialPrecon<CoarseVector> simpleC;
auto RunSolve = [&](int nr)
{
std::cout << GridLogMessage << "**********************************************" << std::endl;
std::cout << GridLogMessage << " V2 THREE-level solve, Nrhs = " << nr << std::endl;
// Device-memory budget BEFORE the solve (2026-08-28: NRHS=6 died in hipMalloc at the
// first fine-smoother history allocation -- reported as an asynchronous "memory
// access fault" unless AMD_SERIALIZE_KERNEL/COPY made it a clean OOM). The outer
// mRHS GCR holds src, sol, r, Az and OuterMmax x (p,q) fine fields PER RHS; the fine
// smoother adds FineSmootherMmax x (p,q) once. The MemoryManager LRU cap
// (--device-mem) must be BELOW what is physically left after the non-LRU allocations
// (comms buffers, dense slab, stencil buffers), or the device fills before anything
// is evicted. Print the estimate, the LRU state, and the device's own free count.
{
uint64_t fieldBytes = (uint64_t)FGrid->lSites()*sizeof(typename LatticeFermionD::scalar_object);
double outerGB = (double)nr*(4 + 2*OuterMmax)*fieldBytes/1.0e9;
double smthGB = (double)(2*FineSmootherMmax + 4)*fieldBytes/1.0e9;
std::cout << GridLogMessage << "Device budget: fine field " << fieldBytes/1.0e6 << " MB; outer GCR history "
<< nr << " x (4 + 2 x " << OuterMmax << ") fields = " << outerGB << " GB; fine smoother history + temps ~ "
<< smthGB << " GB; MemoryManager device LRU " << MemoryManager::DeviceCacheBytes()/1.0e9 << " GB now, cap "
<< MemoryManager::DeviceMaxBytes/1.0e9 << " GB" << std::endl;
// Empty the device LRU: every setup-era Lattice copy (coarse null vectors,
// coarsening temporaries) goes back to host, so the solve's working set
// starts from a clean device and the cap applies to it alone.
MemoryManager::EvictAll();
// ...and release the allocation caches' held blocks (setup-era deviceVector
// scratch that is "free" to the caller but not to hipMalloc).
MemoryManager::DropCache();
MemoryManager::PrintBytes();
acceleratorMem();
}
std::cout << GridLogMessage << "**********************************************" << std::endl;
Coordinate cml({nr,1,clatt[0],clatt[1],clatt[2],clatt[3]});
Coordinate ccml({nr,1,cclatt[0],cclatt[1],cclatt[2],cclatt[3]});
GridCartesian *CMrhs = new GridCartesian(cml, cmsimd,cmmpi);
GridCartesian *CCMrhs = new GridCartesian(ccml,cmsimd,cmmpi);
CoarseOpPV.SetGrid(CMrhs);
CoarseOpL2.SetGrid(CCMrhs);
NonHermitianLinearOperator<CoarseOperator,CoarseVector> LinOpC (CoarseOpPV);
NonHermitianLinearOperator<CoarseCoarseOperator,CoarseCoarseVector> LinOpCC(CoarseOpL2);
MrhsDenseCCSolve<DenseCC_t,CoarseCoarseVector> ccSolve(*DenseCC,nr);
ShiftedLinearOperator<CoarseVector> ShiftedC(CoarseSmootherShift, LinOpC);
if ( PowerIterations > 0 ) {
// Spectral edges the smoother polynomials must respect (see
// scripts/gcr_polynomial.py: |R_m|>1 beyond the edge = amplification).
PowerIteration<CoarseVector> ("CoarseSmootherOp(shift="+std::to_string(CoarseSmootherShift)+")", ShiftedC, CMrhs, PowerIterations);
PowerIteration<CoarseVector> ("CoarseOp(unshifted)", LinOpC, CMrhs, PowerIterations);
PowerIteration<LatticeFermionD>("FineSmootherOp(shift="+std::to_string(FineSmootherShift)+")", ShiftedPVdagM, Ddwf.FermionGrid(), PowerIterations);
}
PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>
CoarseSmootherGCR(0.01,1,ShiftedC,simpleC,CoarseSmootherMmax,CoarseSmootherNstep);
CoarseSmootherGCR.Level(2); CoarseSmootherGCR.Name("Csmoother"); CoarseSmootherGCR.SetZeroGuess(1);
SwitchableSmoother<CoarseVector> CoarseSmootherSlot(CoarseSmootherGCR,"Csmoother GCR");
MrhsCoarseThreeLevelPrec<CoarseVector,CoarseCoarseVector>
L2to3Precon(LinOpC, CoarseSmootherSlot, MrhsProjectorL2, ccSolve,
Coarse5d, CoarseCoarse5d, CCMrhs, nr);
PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>
L2PGCR(CoarseSolverTol, CoarseSolverOrder/16, LinOpC, L2to3Precon, CoarseSolverMmax, 16);
L2PGCR.Level(2); L2PGCR.Name("Couter"); L2PGCR.SetZeroGuess(1);
FineSmoother_t SmootherGCR(0.0,1,ShiftedPVdagM,simple_fine,FineSmootherMmax,FineSmootherOrder);
SmootherGCR.Level(1); SmootherGCR.Name("Fsmoother"); SmootherGCR.SetZeroGuess(1);
SmootherGCR.LogCoefficients(SmootherCoeffLog);
CoarseSmootherGCR.LogCoefficients(SmootherCoeffLog);
SwitchableSmoother<LatticeFermionD> FineSmootherSlot(SmootherGCR,"Fsmoother GCR");
MrhsTwoLevelMG<LatticeFermionD,CoarseVector,SwitchableSmoother<LatticeFermionD> >
ThreeLevelPrecon(PVdagM, FineSmootherSlot, MrhsProjector, L2PGCR, Coarse5d, CMrhs);
ThreeLevelPrecon.SetSloppy = SetFineSloppy;
ThreeLevelPrecon.SloppyComms = FineSloppyComms;
MrhsPGCRNonHermitian<LatticeFermionD>
L1PGCR(OuterTol,1000,PVdagM,ThreeLevelPrecon,OuterMmax,OuterNstep);
L1PGCR.Level(1); L1PGCR.Name("Fouter"); L1PGCR.SetZeroGuess(1);
//////////////////////////////////////////////////////////////////////
// Smoother modes. Objects live for the duration of this RunSolve.
//////////////////////////////////////////////////////////////////////
std::unique_ptr<ChebyshevNonHermitianSmoother<LatticeFermionD> > FineCheb;
std::unique_ptr<ChebyshevNonHermitianSmoother<CoarseVector> > CoarseCheb;
std::unique_ptr<GCRReplaySmoother<LatticeFermionD> > FineReplay;
std::unique_ptr<GCRReplaySmoother<CoarseVector> > CoarseReplay;
GCRCoefficients recF, recC;
{
GCRCoefficients::Select sel = GCRCoefficients::Last;
if ( PolyRecordSelect=="first" ) sel = GCRCoefficients::First;
if ( PolyRecordSelect=="mean" ) sel = GCRCoefficients::Mean;
recF.select = sel; recC.select = sel;
}
if ( FineSmootherMode == "cheb" ) {
FineCheb.reset(new ChebyshevNonHermitianSmoother<LatticeFermionD>(FineChebLo,FineChebHi,FineSmootherOrder,ShiftedPVdagM));
FineCheb->Verbose = PolyVerbose; FineCheb->name = "Fsmoother";
FineSmootherSlot.Set(*FineCheb,"Fsmoother Chebyshev");
}
if ( CoarseSmootherMode == "cheb" ) {
CoarseCheb.reset(new ChebyshevNonHermitianSmoother<CoarseVector>(CoarseChebLo,CoarseChebHi,CoarseSmootherNstep,ShiftedC));
CoarseCheb->Verbose = PolyVerbose; CoarseCheb->name = "Csmoother";
CoarseSmootherSlot.Set(*CoarseCheb,"Csmoother Chebyshev");
}
// Recording window [PolyRecordStart, PolyRecordStart+PolyRecordIters).
// M3 (2026-08-26): the GCR polynomial changes fast over the first outer
// steps (per-call |r|/|r0| 0.0034 -> 0.017 over steps 1-4) and the MEAN
// of those is a poor smoother (replay 0.033-0.049 per call); record a
// settled window instead.
if ( PolyRecordStart == 0 ) {
if ( FineSmootherMode == "replay" ) SmootherGCR.SetCoefficientRecorder(&recF);
if ( CoarseSmootherMode == "replay" ) CoarseSmootherGCR.SetCoefficientRecorder(&recC);
}
// Record -> replay, with optional periodic re-recording ("re-record, not
// fade away": HDCG refreshed its polynomial every 10 steps, tracking the
// evolving spectral content of the residual). Schedule on outer steps:
// [PolyRecordStart, +PolyRecordIters) adaptive GCR, recording
// then replay of the selected recorded call;
// if PolyRefresh>0: every PolyRefresh steps, ONE adaptive recording
// step, then replay of that call.
auto BuildReplays = [&](void){
if ( FineSmootherMode == "replay" ) {
SmootherGCR.SetCoefficientRecorder(nullptr);
recF.Flush(); recF.Report("Fsmoother");
FineReplay.reset(new GCRReplaySmoother<LatticeFermionD>(ShiftedPVdagM,recF));
FineReplay->Verbose = PolyVerbose; FineReplay->name = "Fsmoother";
FineSmootherSlot.Set(*FineReplay,"Fsmoother replay");
SmootherGCR.ReleaseHistory();
}
if ( CoarseSmootherMode == "replay" ) {
CoarseSmootherGCR.SetCoefficientRecorder(nullptr);
recC.Flush(); recC.Report("Csmoother");
CoarseReplay.reset(new GCRReplaySmoother<CoarseVector>(ShiftedC,recC));
CoarseReplay->Verbose = PolyVerbose; CoarseReplay->name = "Csmoother";
CoarseSmootherSlot.Set(*CoarseReplay,"Csmoother replay");
CoarseSmootherGCR.ReleaseHistory();
}
};
auto StartRecording = [&](int step){
if ( FineSmootherMode == "replay" ) { { auto sel=recF.select; recF = GCRCoefficients(); recF.select=sel; } SmootherGCR.SetCoefficientRecorder(&recF); FineSmootherSlot.Set(SmootherGCR,"Fsmoother GCR (recording)"); }
if ( CoarseSmootherMode == "replay" ) { { auto sel=recC.select; recC = GCRCoefficients(); recC.select=sel; } CoarseSmootherGCR.SetCoefficientRecorder(&recC); CoarseSmootherSlot.Set(CoarseSmootherGCR,"Csmoother GCR (recording)"); }
std::cout << GridLogMessage << "Smoother coefficient recording starts at outer step " << step << std::endl;
};
int switchStep = PolyRecordStart + PolyRecordIters;
L1PGCR.OnStep = [&](int step){
if ( FineSmootherMode != "replay" && CoarseSmootherMode != "replay" ) return;
if ( step == PolyRecordStart && PolyRecordStart > 0 ) StartRecording(step);
if ( step == switchStep ) { BuildReplays(); return; }
if ( PolyRefresh > 0 && step > switchStep ) {
int since = step - switchStep;
if ( since % PolyRefresh == 0 ) { StartRecording(step); } // one adaptive, recorded step
if ( since % PolyRefresh == 1 ) { BuildReplays(); } // then replay it
}
};
std::cout << GridLogMessage << "Smoother modes: fine " << FineSmootherMode << " coarse " << CoarseSmootherMode
<< (FineSmootherMode=="replay"||CoarseSmootherMode=="replay" ? " (record outer steps "+std::to_string(PolyRecordStart)+".."+std::to_string(PolyRecordStart+PolyRecordIters)+")" : "")
<< std::endl;
std::vector<LatticeFermionD> src(nr,FGrid), sol(nr,FGrid);
for(int r=0;r<nr;r++){ gaussian(RNG5,src[r]); sol[r]=Zero(); }
GridStopWatch w; w.Start();
L1PGCR(src,sol);
w.Stop();
std::cout << GridLogMessage << "V2 3-level solve Nrhs "<<nr<<" total " << w.Elapsed()
<< " (per RHS: " << w.useconds()/1.0e6/nr << " s)" << std::endl;
// The outer operator is exact by policy; assert the state rather than
// trust it -- a preconditioner that failed to restore would surface here.
SetFineSloppy(0);
{ LatticeFermionD Ax(FGrid); RealD worst=0.0;
for(int r=0;r<nr;r++){ PVdagM.Op(sol[r],Ax); Ax=Ax-src[r];
RealD rn=std::sqrt(norm2(Ax)/norm2(src[r]));
std::cout << GridLogMessage << "FINAL Nrhs "<<nr<<": rhs["<<r<<"] true residual = " << rn << std::endl;
worst=std::max(worst,rn); }
std::cout << GridLogMessage << "FINAL Nrhs "<<nr<<": worst-case residual = " << worst
<< " (exact-halo verification)" << std::endl;
}
// The operators borrow these grids and build a PaddedCell on them, so
// they must let go before the grids are destroyed.
CoarseOpPV.ReleaseGrid();
CoarseOpL2.ReleaseGrid();
delete CMrhs; delete CCMrhs;
};
RunSolve(nrhs);
if ( getenv("SOLVE_SRHS")==nullptr || atoi(getenv("SOLVE_SRHS")) ) RunSolve(1);
std::cout << GridLogMessage << "*** stage three complete: solves done ***" << std::endl;
Grid_finalize();
}
examples/Makefile.am
+9 -3
View File
@@ -2,13 +2,15 @@ SUBDIRS = .
include Make.inc
# Precision variants of the two multigrid drivers. The coarse-sector and
# Compile-time variants of the two multigrid drivers. The coarse-sector and
# dense-inversion precisions are compile-time instantiations, so each is its
# own binary built from the same source; the fine-level precision is a
# run-time parameter and needs no variant. Per-target flags give each
# variant its own object file.
# run-time parameter and needs no variant. The recur variant additionally
# takes the GCR residual norm from the recurrence rather than a reduction.
# Per-target flags give each variant its own object file.
bin_PROGRAMS += Example_pvdagm_multigrid_fp32coarse \
Example_pvdagm_multigrid_fp32coarse_fp32dense \
Example_pvdagm_multigrid_fp32coarse_fp32dense_recur \
Example_hdcg_multigrid_fp32coarse
Example_pvdagm_multigrid_fp32coarse_SOURCES = Example_pvdagm_multigrid.cc
@@ -19,6 +21,10 @@ Example_pvdagm_multigrid_fp32coarse_fp32dense_SOURCES = Example_pvdagm_multigri
Example_pvdagm_multigrid_fp32coarse_fp32dense_CPPFLAGS = -DCOARSE_SINGLE -DGRID_DENSE_INVERSE_SINGLE
Example_pvdagm_multigrid_fp32coarse_fp32dense_LDADD = $(top_builddir)/Grid/libGrid.a
Example_pvdagm_multigrid_fp32coarse_fp32dense_recur_SOURCES = Example_pvdagm_multigrid.cc
Example_pvdagm_multigrid_fp32coarse_fp32dense_recur_CPPFLAGS = -DCOARSE_SINGLE -DGRID_DENSE_INVERSE_SINGLE -DGRID_GCR_RESIDUAL_RECURRENCE
Example_pvdagm_multigrid_fp32coarse_fp32dense_recur_LDADD = $(top_builddir)/Grid/libGrid.a
Example_hdcg_multigrid_fp32coarse_SOURCES = Example_hdcg_multigrid.cc
Example_hdcg_multigrid_fp32coarse_CPPFLAGS = -DCOARSE_SINGLE
Example_hdcg_multigrid_fp32coarse_LDADD = $(top_builddir)/Grid/libGrid.a