From 5372dce0c279c32f6671bc38a8c360578957e0a3 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Wed, 26 Aug 2026 20:15:14 -0400 Subject: [PATCH] Rerecord, not fade away ! --- Grid/algorithms/iterative/GCRCoefficients.h | 58 ++++++++++++---- ...mple_pvdagm_v2_3level_DenseCoarseMatrix.cc | 68 ++++++++++++++++--- systems/Frontier/smoother_modes.job | 49 ++++++++----- tests/debug/Test_poly_smoother.cc | 1 + 4 files changed, 136 insertions(+), 40 deletions(-) diff --git a/Grid/algorithms/iterative/GCRCoefficients.h b/Grid/algorithms/iterative/GCRCoefficients.h index c09d92747..b435ea8f2 100644 --- a/Grid/algorithms/iterative/GCRCoefficients.h +++ b/Grid/algorithms/iterative/GCRCoefficients.h @@ -28,29 +28,57 @@ NAMESPACE_BEGIN(Grid); ////////////////////////////////////////////////////////////////////////////// struct GCRCoefficients { int mmax = 0; - std::vector a_sum; // [k] - std::vector a_n; - std::vector > b_sum; // [k][j] - std::vector > b_n; + // Every recorded call is kept: calls[c] = list of (a_k, [b_kj]) per step. + // A(k)/B(k,j) return the coefficients of the SELECTED call: by default the + // last complete one. Selection "mean" averages coefficients over calls -- + // kept for comparison only: the mean of the coefficients of a nonlinear + // recurrence is not the mean of the polynomials, and in practice (Frontier + // M3, 2026-08-26) it was worse than every individual call. + enum Select { Last=0, First=1, Index=2, Mean=3 }; + Select select = Last; + int index = 0; + typedef std::vector > > Call; + std::vector calls; + Call current; void RecordA(int k, ComplexD a){ - if ( (int)a_sum.size() <= k ) { a_sum.resize(k+1,ComplexD(0.0)); a_n.resize(k+1,0); } - a_sum[k] += a; a_n[k]++; + if ( k==0 && current.size() ) { calls.push_back(current); current.clear(); } + if ( (int)current.size() <= k ) current.resize(k+1); + current[k].first = a; } void RecordB(int k, const std::vector &b){ - if ( (int)b_sum.size() <= k ) { b_sum.resize(k+1); b_n.resize(k+1); } - if ( b_sum[k].size() < b.size() ) { b_sum[k].resize(b.size(),ComplexD(0.0)); b_n[k].resize(b.size(),0); } - for(int j=0;j<(int)b.size();j++){ b_sum[k][j] += b[j]; b_n[k][j]++; } + if ( (int)current.size() <= k ) current.resize(k+1); + current[k].second = b; } - int Steps(void) const { return a_sum.size(); } - int Calls(void) const { return a_n.size() ? a_n[0] : 0; } - ComplexD A(int k) const { return a_sum[k]/(double)a_n[k]; } - int NB(int k) const { return (k<(int)b_sum.size()) ? b_sum[k].size() : 0; } - ComplexD B(int k,int j) const { return b_sum[k][j]/(double)b_n[k][j]; } + void Flush(void){ if ( current.size() ) { calls.push_back(current); current.clear(); } } + int Calls(void) const { return calls.size() + (current.size() ? 1 : 0); } + const Call & Chosen(void) const { + GRID_ASSERT( calls.size() || current.size() ); + if ( calls.empty() ) return current; + if ( select==First ) return calls.front(); + if ( select==Index ) { GRID_ASSERT(index>=0 && index<(int)calls.size()); return calls[index]; } + return calls.back(); + } + int Steps(void) const { return select==Mean ? MeanSteps() : Chosen().size(); } + int NB(int k) const { return select==Mean ? MeanNB(k) : Chosen()[k].second.size(); } + ComplexD A(int k) const { return select==Mean ? MeanA(k) : Chosen()[k].first; } + ComplexD B(int k,int j) const { return select==Mean ? MeanB(k,j) : Chosen()[k].second[j]; } + // mean over calls (comparison only) + int MeanSteps(void) const { int m=0; for(auto &c:calls) m = std::max(m,(int)c.size()); return m; } + int MeanNB(int k) const { int m=0; for(auto &c:calls) if(k<(int)c.size()) m = std::max(m,(int)c[k].second.size()); return m; } + ComplexD MeanA(int k) const { ComplexD s(0.0); int n=0; for(auto &c:calls) if(k<(int)c.size()){ s+=c[k].first; n++; } return s/(double)n; } + ComplexD MeanB(int k,int j) const { ComplexD s(0.0); int n=0; for(auto &c:calls) if(k<(int)c.size() && j<(int)c[k].second.size()){ s+=c[k].second[j]; n++; } return s/(double)n; } void Report(const std::string &name) const { - std::cout << GridLogMessage << "GCRCoefficients " << name << ": " << Calls() << " calls, " << Steps() << " steps, mmax " << mmax << std::endl; + const char *sel[4]={"last","first","index","mean"}; + std::cout << GridLogMessage << "GCRCoefficients " << name << ": " << Calls() << " calls, " << Steps() << " steps, mmax " << mmax + << ", selection " << sel[select] << std::endl; for(int k=0;k1 ) { + RealD lo=1e300, hi=0; for(auto &c:calls) if(k<(int)c.size()){ RealD x=real(c[k].first), y=imag(c[k].first); RealD m=std::sqrt(x*x+y*y); lo=std::min(lo,m); hi=std::max(hi,m); } + std::cout << " |a| over calls [" << lo << "," << hi << "]"; + } std::cout << std::endl; } } diff --git a/examples/Example_pvdagm_v2_3level_DenseCoarseMatrix.cc b/examples/Example_pvdagm_v2_3level_DenseCoarseMatrix.cc index 49d81afb8..a067eed0e 100644 --- a/examples/Example_pvdagm_v2_3level_DenseCoarseMatrix.cc +++ b/examples/Example_pvdagm_v2_3level_DenseCoarseMatrix.cc @@ -103,6 +103,9 @@ int PowerIterations = 0; // >0: power-iterate the smoother operators 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; @@ -153,6 +156,9 @@ void ParseEnvironment(void) 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")); @@ -550,6 +556,17 @@ public: int main (int argc, char ** argv) { + // GRID_MPI_THREAD_MULTIPLE=1: initialise MPI at MPI_THREAD_MULTIPLE before + // Grid_init (Grid asks for SERIALIZED). The SLATE harness does this and its + // copy of the 2D inverse ran at ~2x the production ring rate (62 s vs + // 133-141 s). Second hypothesis behind OMP_NUM_THREADS; test one at a time. + // Pair with MPICH_MAX_THREAD_SAFETY=multiple. + if ( getenv("GRID_MPI_THREAD_MULTIPLE") ) { + int provided = 0; + MPI_Init_thread(&argc, &argv, MPI_THREAD_MULTIPLE, &provided); + std::cout << "GRID_MPI_THREAD_MULTIPLE: requested MPI_THREAD_MULTIPLE, provided " << provided + << (provided==MPI_THREAD_MULTIPLE ? " (MULTIPLE)" : " (NOT multiple)") << std::endl; + } Grid_init(&argc,&argv); ParseEnvironment(); @@ -991,6 +1008,12 @@ int main (int argc, char ** argv) std::unique_ptr > FineReplay; std::unique_ptr > 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(FineChebLo,FineChebHi,FineSmootherOrder,ShiftedPVdagM)); FineCheb->Verbose = PolyVerbose; FineCheb->name = "Fsmoother"; @@ -1001,29 +1024,58 @@ int main (int argc, char ** argv) CoarseCheb->Verbose = PolyVerbose; CoarseCheb->name = "Csmoother"; CoarseSmootherSlot.Set(*CoarseCheb,"Csmoother Chebyshev"); } - if ( FineSmootherMode == "replay" ) SmootherGCR.SetCoefficientRecorder(&recF); - if ( CoarseSmootherMode == "replay" ) CoarseSmootherGCR.SetCoefficientRecorder(&recC); - L1PGCR.OnStep = [&](int step){ - if ( step != PolyRecordIters ) return; + // 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.Report("Fsmoother"); + recF.Flush(); recF.Report("Fsmoother"); FineReplay.reset(new GCRReplaySmoother(ShiftedPVdagM,recF)); FineReplay->Verbose = PolyVerbose; FineReplay->name = "Fsmoother"; FineSmootherSlot.Set(*FineReplay,"Fsmoother replay"); - SmootherGCR.ReleaseHistory(); // memory-neutral swap: the GCR's history goes as the replay's comes + SmootherGCR.ReleaseHistory(); } if ( CoarseSmootherMode == "replay" ) { CoarseSmootherGCR.SetCoefficientRecorder(nullptr); - recC.Report("Csmoother"); + recC.Flush(); recC.Report("Csmoother"); CoarseReplay.reset(new GCRReplaySmoother(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 for "+std::to_string(PolyRecordIters)+" outer steps)" : "") + << (FineSmootherMode=="replay"||CoarseSmootherMode=="replay" ? " (record outer steps "+std::to_string(PolyRecordStart)+".."+std::to_string(PolyRecordStart+PolyRecordIters)+")" : "") << std::endl; std::vector src(nr,FGrid), sol(nr,FGrid); diff --git a/systems/Frontier/smoother_modes.job b/systems/Frontier/smoother_modes.job index befb1dc24..2d56f6421 100644 --- a/systems/Frontier/smoother_modes.job +++ b/systems/Frontier/smoother_modes.job @@ -96,34 +96,49 @@ export PowerIterations=0 export SmootherCoeffLog=0 # frozen-polynomial controls -export PolyRecordIters=4 # outer steps recorded before the switch +export PolyRecordIters=8 # outer steps recorded +export PolyRecordStart=8 # ...starting here: the early-step polynomials are unrepresentative (M3) +export PolyRecordSelect=last # replay ONE recorded call's polynomial (PB: every individual call beats the coefficient mean) +export PolyRefresh=5 # re-record every 5 outer steps: BFM BfmHDCG.C:2243, k%5==1 -> LdopM1MirsPolyRecord, single call, replayed 4 steps +# Inverse ring-rate hypotheses, ONE AT A TIME: (1) OMP_NUM_THREADS=1 (set above); +# (2) if (1) fails, uncomment the two lines below (harness ran 62 s with these). +#export MPICH_MAX_THREAD_SAFETY=multiple +#export GRID_MPI_THREAD_MULTIPLE=1 export PolyVerbose=1 # frozen smoothers print |r_m|/|r_0| per call: separates 'bad polynomial' from 'linear V-cycle stagnates the outer' export FineChebLo=3.0 # harvested |R|<0.1 edge / PowerIteration edge x1.05 export FineChebHi=137.0 export CoarseChebLo=8.0 export CoarseChebHi=47.0 # shift 2.0: edge 43.3 x1.08 -run_mode () { - name=$1; export FineSmootherMode=$2; export CoarseSmootherMode=$3 - echo "----- $name : FineSmootherMode=$FineSmootherMode CoarseSmootherMode=$CoarseSmootherMode -----" - fname=log.modes.$name +# Reference: the banked ADAPTIVE optimum, Fso6 / sm4 / Css2.0 / Nstep2 / svm8 -> +# 28.57 s (Nrhs=1), ~14.9 s/RHS (Nrhs=4). The stationary smoother converged at +# the deliberate overshoot (order 12, fine shift 1.0, coarse Nstep 6); the +# ladder below walks back towards the banked point. A cell wins if it stays +# convergent AND beats 28.57 s. Each cell ~5 min. +run_cell () { + name=$1; export FineSmootherOrder=$2; export FineSmootherShift=$3; export CoarseSmootherNstep=$4 + export FineSmootherMode=$5; export CoarseSmootherMode=$6 + echo "----- $name : Fso=$FineSmootherOrder Fss=$FineSmootherShift Csn=$CoarseSmootherNstep fine=$FineSmootherMode coarse=$CoarseSmootherMode -----" + fname=log.ladder.$name srun -N36 -n288 --kill-on-bad-exit=1 ./select_gpu $root/examples/Example_pvdagm_v2_3level_DenseCoarseMatrix \ --mpi ${MPI_GEOM} --grid $vol $OPTS1 --comms-overlap > $fname 2>&1 - echo " exit $?"; sleep 60 # let a faulted step's GPUs be released before the next srun (M3 after M2 faulted instantly) - echo " $(grep -h 'V2 3-level solve Nrhs' $fname | tr '\n' ' ')" - echo " $(grep -h 'Fouter MrhsPGCR: Converged' $fname | sed 's/.*Converged/Converged/' | tr '\n' ' ')" - echo " $(grep -h 'FINAL Nrhs .: worst' $fname | tr '\n' ' ')" - grep -h "SwitchableSmoother\|GCRCoefficients .*calls" $fname | head -6 + echo " exit $?"; sleep 60 + echo " $(grep -h 'V2 3-level solve Nrhs' $fname | sed 's/.*V2/V2/' | tr '\n' ' ')" + echo " $(grep -h 'Fouter MrhsPGCR: Converged' $fname | sed 's/.*Converged/Converged/' | cut -c1-60 | tr '\n' ' ')" + echo " replay per-call |r|/|r0| (Nrhs=1 solve): $(awk '/THREE-level solve, Nrhs = 1/{s=1} s && /Fsmoother replay \|r\|/{v=$NF; n++; t+=v; if(v>mx)mx=v} END{if(n) printf "mean %.4f max %.4f over %d calls", t/n, mx, n}' $fname)" + grep -h "SCHUR fp64 distributed invert took\|GB/s/rank" $fname | sed 's/^Grid : Message : [0-9.]* s : //' | cut -c1-120 | head -2 } -run_mode M1_gcr_gcr gcr gcr -run_mode M2_replay_replay replay replay -run_mode M3_replay_gcr replay gcr -run_mode M4_cheb_gcr cheb gcr -run_mode M5_cheb_cheb cheb cheb -run_mode M6_gcr_replay gcr replay # coarse frozen only: M2 showed Couter 16 -> 5 steps with it live +# name Fso Fss Csn fine coarse +run_cell L0_overshoot 12 1.0 6 replay gcr # the converged overshoot (M3), now with last-call selection + refresh 5 +run_cell L1_csn2 12 1.0 2 replay gcr # coarse smoother back to the banked 2 steps +run_cell L2_fso8 8 1.0 2 replay gcr +run_cell L3_fss05 8 0.5 2 replay gcr +run_cell L4_banked 6 0.5 2 replay gcr # nearest to the banked adaptive point +run_cell L5_both 8 0.5 2 replay replay # coarse frozen too, at the best-looking fine point echo "=========================================================" echo "summary" -grep -h "V2 3-level solve Nrhs 1" log.modes.* | sed 's/.*V2/V2/' +for f in log.ladder.*; do echo "$f: $(grep -h "V2 3-level solve Nrhs 1" $f | sed "s/.*V2/V2/") outer $(grep -h "Fouter MrhsPGCR: Converged" $f | tail -1 | grep -oE "iteration [0-9]+")"; done +echo "reference (adaptive, banked): 28.57 s Nrhs=1, 14.9 s/RHS Nrhs=4" echo "=========================================================" diff --git a/tests/debug/Test_poly_smoother.cc b/tests/debug/Test_poly_smoother.cc index b3b7b2dae..36d89ba9f 100644 --- a/tests/debug/Test_poly_smoother.cc +++ b/tests/debug/Test_poly_smoother.cc @@ -104,6 +104,7 @@ int main(int argc, char **argv) LatticeFermionD src(UGrid), x(UGrid); for(int c=0;c