Both changes to improve the buffer persistence in teh dense distributed inverse and to add verbosity to polynomial smoothers

Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h

@@ -85,6 +85,9 @@ public:
   // the iteration.
   GCRCoefficients *Recorder = nullptr;
   void SetCoefficientRecorder(GCRCoefficients *r) { Recorder = r; if(r) r->mmax = mmax; };
+  // Free the persistent history (e.g. when this solver is replaced by a
+  // replayed polynomial): re-made on the next call if ever needed again.
+  void ReleaseHistory(void) { q.clear(); p.clear(); qq.clear(); hist_grid = nullptr; };
 
   PrecGeneralisedConjugateResidualNonHermitian(RealD tol,Integer maxit,LinearOperatorBase<Field> &_Linop,LinearFunction<Field> &Prec,int _mmax,int _nstep) : 
     Tolerance(tol), 

Grid/algorithms/multigrid/BlockCyclic.h

@@ -62,13 +62,19 @@ public:
   // For fixed P the per-rank SUMMA volume  N^2 (1/Pr + 1/Pc)  is minimised
   // at the most square grid.  P=288 -> 16 x 18.
   ///////////////////////////////////////////////////////////////////////////
+  // Nearest-to-square factorisation with Pr >= Pc.  The orientation matters
+  // for SUMMA: B panels travel along process columns in slots of
+  // S = ceil(Pc/Pr) panels, padded; with Pr >= Pc, S == 1 and there is no
+  // padding.  (288 as 16x18 had S=2 -- half the B-ring bytes were zeros;
+  // 18x16 moves ~30% fewer bytes for the same inverse.)
   static void ChooseProcessGrid(int P, int &Pr, int &Pc)
   {
     GRID_ASSERT(P >= 1);
-    Pr = 1;
+    int r = 1;
-    for(int r=1; (int64_t)r*r <= (int64_t)P; r++)
+    for(int f=1; (int64_t)f*f <= (int64_t)P; f++)
-      if ( P % r == 0 ) Pr = r;
+      if ( P % f == 0 ) r = f;
-    Pc = P / Pr;
+    Pc = r;
+    Pr = P / r;
   }
 
   ///////////////////////////////////////////////////////////////////////////

Grid/algorithms/multigrid/BlockCyclicSumma.h

@@ -132,6 +132,15 @@ public:
   // row, ring B along the process column, then the local GEMMs.  The rings
   // are synchronous SendToRecvFrom, so tRing is time the GPU is idle unless
   // a future version overlaps them with the GEMMs.
+  // PERSISTENT ring buffers: allocated once (grow-only) and reused by every
+  // Multiply, so the device addresses handed to MPI never change.  Fresh
+  // per-call buffers rotated through the caching allocator's blocks, and a
+  // GPU-direct RDMA registration cache keyed on address then re-registers
+  // per message: measured 1.26 GB/s/rank on 13 MB ring messages (production,
+  // GRID_ALLOC_NCACHE_LARGE=64) against 62 s total for the same inverse when
+  // hipMalloc returned a stable address.
+  deviceVector<ComplexD> Abuf;
+  deviceVector<ComplexD> Bbuf;
   double   tAlloc=0, tPack=0, tRingA=0, tRingB=0, tGemm=0;
   uint64_t bytesRing=0, nRingMsg=0, nMultiply=0, nGemm=0;
   void ResetTelemetry(void){ tAlloc=tPack=tRingA=tRingB=tGemm=0; bytesRing=nRingMsg=nMultiply=nGemm=0; }
@@ -199,8 +208,8 @@ public:
     const uint64_t slotB  = (uint64_t)S*slotB1;
     nMultiply++;
     tAlloc -= usecond();
-    deviceVector<ComplexD> Abuf( slotA*Pc ? slotA*Pc : 1 );
+    if ( Abuf.size() < std::max<uint64_t>(slotA*Pc,1) ) Abuf.resize( std::max<uint64_t>(slotA*Pc,1) );
-    deviceVector<ComplexD> Bbuf( slotB*Pr ? slotB*Pr : 1 );
+    if ( Bbuf.size() < std::max<uint64_t>(slotB*Pr,1) ) Bbuf.resize( std::max<uint64_t>(slotB*Pr,1) );
     tAlloc += usecond();
 
     deviceVector<ComplexD *> ap(1), bp(1), cp(1);

Grid/algorithms/multigrid/Smoothers.h

@@ -91,6 +91,8 @@ public:
   LinearOperatorBase<Field> &Linop;
   RealD lo, hi; int order;
   std::vector<RealD> Coeffs;
+  int Verbose = 0;
+  std::string name = "cheb";
   ChebyshevNonHermitianSmoother(RealD _lo,RealD _hi,int _order,LinearOperatorBase<Field> &Op)
     : Linop(Op), lo(_lo), hi(_hi), order(_order)
   {
@@ -126,6 +128,7 @@ public:
       if ( Coeffs[n] != 0.0 ) axpy(out,Coeffs[n],*Tnp,out);
       Field *swizzle=Tnm; Tnm=Tn; Tn=Tnp; Tnp=swizzle;
     }
+    if ( Verbose ) { Linop.Op(out,y); y = y - in; std::cout << GridLogMessage << " " << name << " cheb |r|/|r0| = " << std::sqrt(norm2(y)/norm2(in)) << std::endl; }
   }
 };
 
@@ -205,16 +208,27 @@ public:
   std::vector<std::vector<ComplexD> > b;
   GridBase *hist_grid = nullptr;
   std::vector<Field> p;
+  int Verbose = 0;                       // 1: print |r_m|/|r_0| per call (one extra reduction)
+  std::string name = "replay";
   GCRReplaySmoother(LinearOperatorBase<Field> &Op, const GCRCoefficients &c) : Linop(Op)
   {
     mmax  = c.mmax;  GRID_ASSERT(mmax>=1);
     nstep = c.Steps(); GRID_ASSERT(nstep>=1);
     a.resize(nstep); b.resize(nstep);
+    RealD amax=0.0, bmax=0.0;
     for(int k=0;k<nstep;k++){
       a[k] = c.A(k);
+      GRID_ASSERT( std::isfinite(real(a[k])) && std::isfinite(imag(a[k])) );
+      amax = std::max(amax, std::abs(a[k]));
       b[k].resize(c.NB(k));
-      for(int j=0;j<c.NB(k);j++) b[k][j] = c.B(k,j);
+      GRID_ASSERT( c.NB(k) <= mmax-1 );
+      for(int j=0;j<c.NB(k);j++){
+        b[k][j] = c.B(k,j);
+        GRID_ASSERT( std::isfinite(real(b[k][j])) && std::isfinite(imag(b[k][j])) );
+        bmax = std::max(bmax, std::abs(b[k][j]));
+      }
     }
+    std::cout << GridLogMessage << " GCRReplaySmoother: max|a| " << amax << " max|b| " << bmax << std::endl;
     std::cout << GridLogMessage << " GCRReplaySmoother: " << nstep << " steps, mmax " << mmax
               << ", from " << c.Calls() << " recorded calls" << std::endl;
   }
@@ -231,11 +245,15 @@ public:
     r = src;
     psi = Zero();
     p[0] = r;
+    RealD r0 = Verbose ? norm2(src) : 0.0;
     for(int k=0;k<nstep;k++){
       int kp=k+1, peri_k=k%mmax, peri_kp=kp%mmax;
       Linop.Op(p[peri_k],q);                 // q_k = A p_k
       axpy(psi,  a[k], p[peri_k], psi);
-      if ( k==nstep-1 ) break;
+      if ( k==nstep-1 ) {
+        if ( Verbose ) { axpy(r,-a[k],q,r); std::cout << GridLogMessage << " " << name << " replay |r|/|r0| = " << std::sqrt(norm2(r)/r0) << std::endl; }
+        break;
+      }
       axpy(r,   -a[k], q,         r);
       p[peri_kp] = r;
       for(int j=0;j<(int)b[k].size();j++){

examples/Example_pvdagm_v2_3level_DenseCoarseMatrix.cc

@@ -103,6 +103,7 @@ int   PowerIterations      = 0;     // >0: power-iterate the smoother operators
 std::string FineSmootherMode   = "gcr";
 std::string CoarseSmootherMode = "gcr";
 int   PolyRecordIters      = 4;
+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;
@@ -152,6 +153,7 @@ void ParseEnvironment(void)
   if(getenv("FineSmootherMode"))   FineSmootherMode   = getenv("FineSmootherMode");
   if(getenv("CoarseSmootherMode")) CoarseSmootherMode = getenv("CoarseSmootherMode");
   if(getenv("PolyRecordIters"))    PolyRecordIters    = atoi(getenv("PolyRecordIters"));
+  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"));
@@ -991,10 +993,12 @@ int main (int argc, char ** argv)
     GCRCoefficients recF, recC;
     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");
     }
     if ( FineSmootherMode   == "replay" ) SmootherGCR.SetCoefficientRecorder(&recF);
@@ -1005,13 +1009,17 @@ int main (int argc, char ** argv)
         SmootherGCR.SetCoefficientRecorder(nullptr);
         recF.Report("Fsmoother");
         FineReplay.reset(new GCRReplaySmoother<LatticeFermionD>(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
       }
       if ( CoarseSmootherMode == "replay" ) {
         CoarseSmootherGCR.SetCoefficientRecorder(nullptr);
         recC.Report("Csmoother");
         CoarseReplay.reset(new GCRReplaySmoother<CoarseVector>(ShiftedC,recC));
+        CoarseReplay->Verbose = PolyVerbose; CoarseReplay->name = "Csmoother";
         CoarseSmootherSlot.Set(*CoarseReplay,"Csmoother replay");
+        CoarseSmootherGCR.ReleaseHistory();
       }
     };
     std::cout << GridLogMessage << "Smoother modes: fine " << FineSmootherMode << "  coarse " << CoarseSmootherMode

systems/Frontier/schur2d_ladder.job

@@ -80,7 +80,7 @@ S2D_N=13824 srun -N1 -n8 ./select_gpu $root/tests/debug/Test_schur2d_scale \
 
 ##############################################################################
 echo "========================================================="
-echo "F3: scale rehearsal, 36 nodes, N=138240 (nb=480, grid 16x18)"
+echo "F3: scale rehearsal, 36 nodes, N=138240 (nb=480, grid 18x16)"
 echo "     THE number: invert phase vs the 1D baseline 328 s"
 echo "========================================================="
 ##############################################################################

systems/Frontier/slate_vs_schur2d.job

@@ -20,7 +20,7 @@
 #
 #   S1 : 1 node,  8 ranks, N=4096    -- API/device/thread-level shakeout
 #   S2 : 1 node,  8 ranks, N=13824   -- one-node production-shape rehearsal
-#   S3 : 36 nodes, 288 ranks, N=138240 (nb=480, grid 16x18) -- THE comparison
+#   S3 : 36 nodes, 288 ranks, N=138240 (nb=480, grid 18x16) -- THE comparison
 #
 # S3 is gated on S1 passing (rc=0 and both certificates printed): a broken
 # SLATE leg should not burn the 36-node allocation.
@@ -100,7 +100,7 @@ S2D_N=13824 srun -N1 -n8 ./select_gpu $BIN --mpi 1.1.2.4 --grid 16.16.16.16 $OPT
 
 ##############################################################################
 echo "========================================================="
-echo "S3: THE comparison, 36 nodes, 288 ranks, N=138240 (nb=480, grid 16x18)"
+echo "S3: THE comparison, 36 nodes, 288 ranks, N=138240 (nb=480, grid 18x16)"
 echo "========================================================="
 ##############################################################################
 if [ "$S1RC" -eq 0 ] && [ "$S1CERT" -ge 2 ]

systems/Frontier/smoother_modes.job

@@ -23,6 +23,7 @@
 #                        heavy one at Nrhs=1)
 #   M4 cheb   / gcr      fine Chebyshev [FineChebLo,FineChebHi] order Fso
 #   M5 cheb   / cheb
+#   M6 gcr    / replay   coarse frozen only
 #
 # Laptop 8^4 findings (hot config, Ls=4, NBASIS=8): replay/replay converges
 # (28 vs 23 outer); cheb on the FINE level diverges there while cheb on the
@@ -96,6 +97,7 @@ export SmootherCoeffLog=0
 
 # frozen-polynomial controls
 export PolyRecordIters=4           # outer steps recorded before the switch
+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
@@ -105,8 +107,9 @@ run_mode () {
   name=$1; export FineSmootherMode=$2; export CoarseSmootherMode=$3
   echo "----- $name : FineSmootherMode=$FineSmootherMode CoarseSmootherMode=$CoarseSmootherMode -----"
   fname=log.modes.$name
-  srun -N36 -n288 ./select_gpu $root/examples/Example_pvdagm_v2_3level_DenseCoarseMatrix \
+  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' ' ')"
@@ -118,6 +121,7 @@ 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
 
 echo "========================================================="
 echo "summary"

tests/debug/Test_blockcyclic.cc

@@ -38,7 +38,7 @@ Author: Peter Boyle <pboyle@bnl.gov>
 //        equal mloc*nloc; LocalOffset is a bijection onto [0,mloc*nloc).
 //   T5 : block contiguity: within any owned global block, consecutive
 //        global rows are consecutive local rows (what SUMMA panels rely on).
-//   T6 : ChooseProcessGrid: exact factorisation, Pr<=Pc, most-square.
+//   T6 : ChooseProcessGrid: exact factorisation, Pr>=Pc, most-square.
 //////////////////////////////////////////////////////////////////////////////
 
 #include <Grid/Grid.h>
@@ -227,15 +227,15 @@ int main(int argc, char **argv)
       int Pr,Pc;
       BlockCyclicLayout::ChooseProcessGrid(P,Pr,Pc);
       if ( Pr*Pc != P ) ok = false;
-      if ( Pr > Pc )    ok = false;
+      if ( Pr < Pc )    ok = false;                 // Pr >= Pc (SUMMA B panels unpadded)
-      // most-square: no divisor r with Pr < r <= sqrt(P)
+      // most-square: no divisor r with Pc < r <= sqrt(P)
-      for(int r=Pr+1; (int64_t)r*r <= (int64_t)P; r++)
+      for(int r=Pc+1; (int64_t)r*r <= (int64_t)P; r++)
         if ( P % r == 0 ) ok = false;
     }
     int Pr,Pc;
     BlockCyclicLayout::ChooseProcessGrid(288,Pr,Pc);
-    if ( !(Pr==16 && Pc==18) ) ok = false;
+    if ( !(Pr==18 && Pc==16) ) ok = false;   // Pr >= Pc: SUMMA B panels unpadded (S = ceil(Pc/Pr) = 1)
-    Report("T6  ChooseProcessGrid exact, Pr<=Pc, most-square (288 -> 16x18)", ok);
+    Report("T6  ChooseProcessGrid exact, Pr>=Pc, most-square (288 -> 18x16)", ok);
   }
 
   ////////////////////////////////////////////////////////////////////////

tests/debug/Test_poly_smoother.cc

@@ -0,0 +1,200 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/debug/Test_poly_smoother.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.
+
+    See the full license in the file "LICENSE" in the top level distribution
+    directory
+*************************************************************************************/
+/*  END LEGAL */
+
+//////////////////////////////////////////////////////////////////////////////
+// Fixed-polynomial smoothers (Smoothers.h) against the adaptive GCR they
+// replace, on a shifted Wilson operator (non-Hermitian, spectrum in the right
+// half plane):
+//
+//   record : GCR(mmax=2, 8 steps) with a coefficient recorder on 16 sources
+//   T1     : GCRReplaySmoother on a FRESH source reduces the residual to
+//            within a factor 2 of the live GCR on the same source
+//   T2     : replay and GCR solutions agree to the coefficient spread (<10%)
+//   T3     : ChebyshevNonHermitianSmoother with Op=HermOp reproduces the legacy
+//            HermOp ChebyshevSmoother (real spectrum); on the complex-spectrum
+//            Wilson op its (poor) reduction is printed for information only
+//   T4     : the historical HermOp ChebyshevSmoother compiles with both
+//            constructor signatures and reduces the MdagM residual
+//   T5     : replay is bitwise repeatable (no reductions => no reordering)
+//
+//   mpirun -n 1 ./Test_poly_smoother --grid 8.8.8.8 --mpi 1.1.1.1
+//////////////////////////////////////////////////////////////////////////////
+
+#include <Grid/Grid.h>
+#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h>
+
+using namespace Grid;
+
+static int failures = 0;
+static void Report(const std::string &name, bool pass, const std::string &detail="")
+{
+  std::cout << GridLogMessage << "  " << name << (pass ? "  PASS" : "  ** FAIL **");
+  if ( detail.size() ) std::cout << "   " << detail;
+  std::cout << std::endl;
+  if ( !pass ) failures++;
+}
+
+template<class Field>
+class ShiftedOp : public LinearOperatorBase<Field> {
+  LinearOperatorBase<Field> &_Op; RealD shift;
+public:
+  ShiftedOp(RealD s, LinearOperatorBase<Field> &Op) : _Op(Op), shift(s) {}
+  void OpDiag  (const Field &in, Field &out) { GRID_ASSERT(0); }
+  void OpDir   (const Field &in, Field &out,int dir,int disp) { GRID_ASSERT(0); }
+  void OpDirAll(const Field &in, std::vector<Field> &out) { GRID_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){ GRID_ASSERT(0); }
+  void HermOp  (const Field &in, Field &out) { Field tmp(in.Grid()); Op(in,tmp); AdjOp(tmp,out); }
+};
+
+template<class Field>
+RealD Residual(LinearOperatorBase<Field> &Op, const Field &src, const Field &x)
+{
+  Field r(src.Grid()); Op.Op(x,r); r = r - src;
+  return std::sqrt(norm2(r)/norm2(src));
+}
+
+int main(int argc, char **argv)
+{
+  Grid_init(&argc, &argv);
+
+  GridCartesian         *UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(),
+                                       GridDefaultSimd(Nd, vComplexD::Nsimd()), GridDefaultMpi());
+  GridRedBlackCartesian *UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+
+  std::vector<int> seeds({1,2,3,4});
+  GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds);
+
+  LatticeGaugeFieldD Umu(UGrid);
+  SU<Nc>::HotConfiguration(RNG4, Umu);
+
+  RealD mass = 0.5, shift = 0.1;
+  WilsonFermionD Dw(Umu, *UGrid, *UrbGrid, mass);
+  NonHermitianLinearOperator<WilsonFermionD, LatticeFermionD> Op(Dw);
+  ShiftedOp<LatticeFermionD> SOp(shift, Op);
+  TrivialPrecon<LatticeFermionD> simple;
+
+  const int mmax = 2, nstep = 8, ncal = 16;
+
+  //////////////////////////////////////////////////////////////////////
+  // record
+  //////////////////////////////////////////////////////////////////////
+  PrecGeneralisedConjugateResidualNonHermitian<LatticeFermionD> GCR(0.0, 1, SOp, simple, mmax, nstep);
+  GCR.SetZeroGuess(1); GCR.Name("smoother");
+  GCRCoefficients rec;
+  GCR.SetCoefficientRecorder(&rec);
+  LatticeFermionD src(UGrid), x(UGrid);
+  for(int c=0;c<ncal;c++){ gaussian(RNG4,src); x = Zero(); GCR(src,x); }
+  GCR.SetCoefficientRecorder(nullptr);
+  rec.Report("smoother");
+  Report("record: steps and calls", rec.Steps()==nstep && rec.Calls()==ncal,
+         std::to_string(rec.Steps())+" steps, "+std::to_string(rec.Calls())+" calls");
+
+  //////////////////////////////////////////////////////////////////////
+  // fresh source: live GCR vs replay
+  //////////////////////////////////////////////////////////////////////
+  gaussian(RNG4,src);
+  LatticeFermionD xg(UGrid), xr(UGrid), xr2(UGrid), xc(UGrid), d(UGrid);
+  xg = Zero(); GCR(src,xg);
+  RealD rg = Residual(SOp,src,xg);
+
+  GCRReplaySmoother<LatticeFermionD> Replay(SOp, rec);
+  Replay(src,xr);
+  RealD rr = Residual(SOp,src,xr);
+  Report("T1 replay residual within 2x of live GCR", rr < 2.0*rg,
+         "GCR |r|/|r0| = "+std::to_string(rg)+"  replay "+std::to_string(rr));
+
+  d = xr - xg;
+  RealD rel = std::sqrt(norm2(d)/norm2(xg));
+  Report("T2 replay solution vs GCR solution", rel < 0.1, "rel "+std::to_string(rel));
+
+  //////////////////////////////////////////////////////////////////////
+  // Chebyshev 1/x on [lo,hi], hi from a power iteration on SOp
+  //////////////////////////////////////////////////////////////////////
+  RealD hi;
+  {
+    LatticeFermionD v(UGrid), Av(UGrid); gaussian(RNG4,v);
+    RealD n = std::sqrt(norm2(v)); v = v*(1.0/n);
+    for(int i=0;i<60;i++){ SOp.Op(v,Av); hi = std::sqrt(norm2(Av)); v = Av*(1.0/hi); }
+  }
+  RealD lo = 0.5;
+  std::cout << GridLogMessage << "power iteration |lambda_max| ~ " << hi << "  Chebyshev range [" << lo << "," << 1.05*hi << "]" << std::endl;
+  ChebyshevNonHermitianSmoother<LatticeFermionD> Cheb(lo, 1.05*hi, nstep, SOp);
+  Cheb(src,xc);
+  RealD rc = Residual(SOp,src,xc);
+  // Not gated: Wilson's spectrum has O(1) imaginary parts and a real-interval
+  // Chebyshev fit to 1/x degrades exponentially off the axis (Bernstein
+  // ellipse).  Printed as the reminder of what a non-real spectrum does.
+  std::cout << GridLogMessage << "  (info) ChebyshevNonHermitian on the COMPLEX-spectrum Wilson op: |r|/|r0| = " << rc
+            << "   (GCR " << rg << ")  -- expected poor; PVdagM smoother ops have real coefficients" << std::endl;
+
+  //////////////////////////////////////////////////////////////////////
+  // legacy HermOp Chebyshev smoother, both constructor forms, on MdagM
+  //////////////////////////////////////////////////////////////////////
+  {
+    MdagMLinearOperator<WilsonFermionD, LatticeFermionD> HermOp(Dw);
+    LatticeFermionD hsrc(UGrid), hx(UGrid), hr(UGrid); gaussian(RNG4,hsrc);
+    RealD hhi = 0.0;
+    { LatticeFermionD v(UGrid), Av(UGrid); gaussian(RNG4,v); RealD n=std::sqrt(norm2(v)); v=v*(1.0/n);
+      for(int i=0;i<40;i++){ HermOp.HermOp(v,Av); hhi=std::sqrt(norm2(Av)); v=Av*(1.0/hhi); } }
+    ChebyshevSmoother<LatticeFermionD>                 S4(0.5, 1.05*hhi, 12, HermOp);
+    ChebyshevSmoother<LatticeFermionD, WilsonFermionD> S5(0.5, 1.05*hhi, 12, HermOp, Dw);   // historical 5-arg form
+    S4(hsrc,hx); HermOp.HermOp(hx,hr); hr = hr - hsrc;
+    RealD r4 = std::sqrt(norm2(hr)/norm2(hsrc));
+    S5(hsrc,hx); HermOp.HermOp(hx,hr); hr = hr - hsrc;
+    RealD r5 = std::sqrt(norm2(hr)/norm2(hsrc));
+    Report("T4 legacy ChebyshevSmoother (4- and 5-arg) reduces MdagM residual", r4 < 0.5 && r5 == r4,
+           "|r|/|r0| = "+std::to_string(r4)+" / "+std::to_string(r5));
+
+    // T3: the Op()-based Clenshaw on a REAL-spectrum operator must reproduce
+    // the legacy HermOp smoother: same coefficients, same recurrence.
+    struct HermAsOp : public LinearOperatorBase<LatticeFermionD> {
+      LinearOperatorBase<LatticeFermionD> &H;
+      HermAsOp(LinearOperatorBase<LatticeFermionD> &h) : H(h) {}
+      void OpDiag  (const LatticeFermionD &in, LatticeFermionD &out) { GRID_ASSERT(0); }
+      void OpDir   (const LatticeFermionD &in, LatticeFermionD &out,int dir,int disp) { GRID_ASSERT(0); }
+      void OpDirAll(const LatticeFermionD &in, std::vector<LatticeFermionD> &out) { GRID_ASSERT(0); }
+      void Op      (const LatticeFermionD &in, LatticeFermionD &out) { H.HermOp(in,out); }
+      void AdjOp   (const LatticeFermionD &in, LatticeFermionD &out) { H.HermOp(in,out); }
+      void HermOpAndNorm(const LatticeFermionD &in, LatticeFermionD &out,RealD &n1,RealD &n2){ GRID_ASSERT(0); }
+      void HermOp  (const LatticeFermionD &in, LatticeFermionD &out) { H.HermOp(in,out); }
+    } HOp(HermOp);
+    ChebyshevNonHermitianSmoother<LatticeFermionD> C3(0.5, 1.05*hhi, 12, HOp);
+    LatticeFermionD hx3(UGrid), dd(UGrid);
+    C3(hsrc,hx3); HermOp.HermOp(hx3,hr); hr = hr - hsrc;
+    RealD r3 = std::sqrt(norm2(hr)/norm2(hsrc));
+    S4(hsrc,hx); dd = hx - hx3;
+    RealD reldiff = std::sqrt(norm2(dd)/norm2(hx));
+    Report("T3 ChebyshevNonHermitian(Op=HermOp) == legacy ChebyshevSmoother", r3 < 0.5 && reldiff < 1.0e-12,
+           "|r|/|r0| = "+std::to_string(r3)+"  rel diff of solutions "+std::to_string(reldiff));
+  }
+
+  //////////////////////////////////////////////////////////////////////
+  // determinism
+  //////////////////////////////////////////////////////////////////////
+  Replay(src,xr2);
+  d = xr - xr2;
+  Report("T5 replay bitwise repeatable", norm2(d)==0.0);
+
+  std::cout << GridLogMessage << (failures ? "Test_poly_smoother: FAILURES" : "Test_poly_smoother: ALL PASS") << std::endl;
+  Grid_finalize();
+  return failures ? 1 : 0;
+}