FOrmat warn avoidance

Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h

@@ -57,6 +57,11 @@ public:
   GridStopWatch LinalgTimer;
   std::string name;
   int ZeroGuess  = 0;  // caller contract: guess is always zero => first-cycle r0 = src, skip the apply
+  // persistent GCR history (see GCRnStep)
+  GridBase           *hist_grid = nullptr;
+  std::vector<Field>  q;
+  std::vector<Field>  p;
+  std::vector<RealD>  qq;
   int FirstCycle = 0;
 
   LinearFunction<Field>     &Preconditioner;
@@ -67,6 +72,13 @@ public:
   void Level(int n) { Name("Level " + std::to_string(n)); level = n; }
 
   void SetZeroGuess(int z) { ZeroGuess = z; };
+  // Coefficient logging: one line per step with the step length a_k and the
+  // orthogonalisation coefficients b_{k,j}.  These are the data from which a
+  // FIXED polynomial smoother can be harvested: if they are stable from call
+  // to call, the adaptive GCR can be replaced by a stationary p(A) with the
+  // same applies and no reductions.  Off by default; boss rank prints.
+  int  LogCoeffs = 0;
+  void LogCoefficients(int l) { LogCoeffs = l; };
 
   PrecGeneralisedConjugateResidualNonHermitian(RealD tol,Integer maxit,LinearOperatorBase<Field> &_Linop,LinearFunction<Field> &Prec,int _mmax,int _nstep) : 
     Tolerance(tol), 
@@ -145,9 +157,20 @@ public:
     ////////////////////////////////
     // history for flexible orthog
     ////////////////////////////////
-    std::vector<Field> q(mmax,grid);
-    std::vector<Field> p(mmax,grid);
-    std::vector<RealD> qq(mmax);
+    // History arrays are PERSISTENT across calls (allocated once per grid,
+    // re-made only if the grid or mmax changes).  The per-call form
+    // std::vector<Field>(mmax,grid) built a temporary and copy-constructed it
+    // mmax times on every restart cycle -- measured ~5 ms per fine-smoother
+    // call.  Safe: every entry is written before it is read within a cycle
+    // (q[kp],p[kp] assigned before the northog loop can reach them), so no
+    // stale content is ever consumed.
+    if ( hist_grid != grid || (int)q.size() != mmax ) {
+      q.clear(); p.clear();
+      q.reserve(mmax); p.reserve(mmax);
+      for(int i=0;i<mmax;i++){ q.emplace_back(grid); p.emplace_back(grid); }
+      qq.assign(mmax, 0.0);
+      hist_grid = grid;
+    }
       
     GCRLogLevel<< "PGCR nStep("<<nstep<<")"<<std::endl;
 
@@ -198,7 +221,7 @@ public:
     GCRLogLevel<< "PGCR true residual "<< sqrt(cp/SSQ)     <<std::endl;
 
     for(int k=0;k<nstep;k++){
-
+      GRID_TRACE("PGCR_step");
       steps++;
 
       int kp     = k+1;
@@ -213,6 +236,10 @@ public:
 
       cp = axpy_norm(r,-a,q[peri_k],r);
       LinalgTimer.Stop();
+      if ( LogCoeffs ) {
+        GCRLogLevel<<"coeff["<<k<<"] a=("<<real(a)<<","<<imag(a)<<")"
+                   <<" |r|/|r0|="<<sqrt(cp/SSQ)<<std::endl;
+      }
 
       GCRLogLevel<< "PGCR step["<<steps<<"]  resid " << sqrt(cp/SSQ)<<std::endl;
 
@@ -236,6 +263,7 @@ public:
       p[peri_kp]=z;
 
       int northog = ((kp)>(mmax-1))?(mmax-1):(kp);  // if more than mmax done, we orthog all mmax history.
+      std::ostringstream bs;
       for(int back=0;back<northog;back++){
 
 	int peri_back=(k-back)%mmax;   	  GRID_ASSERT((k-back)>=0);
@@ -243,8 +271,12 @@ public:
 	b=-real(innerProduct(q[peri_back],Az))/qq[peri_back];
 	p[peri_kp]=p[peri_kp]+b*p[peri_back];
 	q[peri_kp]=q[peri_kp]+b*q[peri_back];
+        if ( LogCoeffs ) bs<<" b["<<back<<"]="<<b;
 
       }
+      if ( LogCoeffs && northog ) {
+        GCRLogLevel<<"coeff["<<k<<"]"<<bs.str()<<std::endl;
+      }
       qq[peri_kp]=norm2(q[peri_kp]); // could use axpy_norm
       LinalgTimer.Stop();
     }

Grid/algorithms/multigrid/RecursiveSchurInverse.h

@@ -120,7 +120,9 @@ public:
   uint64_t              bytesAllreduce;
   uint64_t              nAllreduce;    // panel collectives
   uint64_t              nGatherGemm;   // GatherGemm calls
-  uint64_t              nGather;       // AllGatherV collectives (debug gate)
+  uint64_t              nGather;       // gather-path collectives (debug gate)
+  uint64_t              nGatherV;      // ... of which MPI_Allgatherv
+  uint64_t              nBcast;        // ... of which Bcast-assembled
 
   // Persistent grow-only device panel; assembly and collectives are
   // device-resident.  Device builds require GPU-aware MPI.
@@ -368,7 +370,10 @@ public:
       uint64_t panelWords = (uint64_t)k*nchunk;
       uint64_t panelBytesThis = panelWords*sizeof(ComplexD);
 
-      int gatherThis = useGather && ( (int64_t)panelBytesThis >= gatherMinBytes );
+      // The MODE (0/1/2), not a boolean: `useGather && cond` collapses 2 to 1
+      // and silently dispatched DENSE_GATHER=2 onto the known-broken
+      // AllGatherV path (Frontier hang, 2026-08-24).  Identical on all ranks.
+      int gatherThis = ( (int64_t)panelBytesThis >= gatherMinBytes ) ? useGather : 0;
 
       // Arrival skew is charged to the barrier, so that tAllreduce measures
       // transfer alone.  Diagnostic only; off by default.
@@ -417,6 +422,7 @@ public:
         }
         if ( gatherThis == 1 )
         {
+          nGatherV++;
           void *send = owner ? (void *)B.ColumnWindow(colB+j0) : (void *)&dRecvBuf[0];
           grid->AllGatherV(send, counts[me],
                            (void *)&dRecvBuf[0], counts, displs, sizeof(ComplexD));
@@ -431,6 +437,7 @@ public:
           // Every rank issues all C broadcasts in the same order, so
           // collective matching stays positional and safe.
           ////////////////////////////////////////////////////////////////
+          nBcast++;
           if ( owner )
           {
             int64_t off_me = rowStart[me] - rowStart[rB0];
@@ -777,6 +784,8 @@ public:
     nAllreduce     = 0;
     nGatherGemm    = 0;
     nGather        = 0;
+    nGatherV       = 0;
+    nBcast         = 0;
 
     SchurNode(0, P, 0, N, Arows);
 
@@ -820,7 +829,8 @@ public:
               << "  leaf "      << tLeaf/1.0e6
               << std::endl;
     std::cout << GridLogMessage << "Schur comms:"
-              << ( useGather ? "  [AllGatherV]" : "  [zero-fill+GlobalSum]" )
+              << "  [transports run: allreduce " << (nAllreduce - nGatherV - nBcast)
+              << " allgatherv " << nGatherV << " bcast " << nBcast << "]"
               << ( barrierProbe ? " [barrier probe on]" : "" )
               << "  GatherGemm calls " << nGatherGemm
               << "  panel allreduces " << nAllreduce

Grid/util/FlightRecorder.cc

@@ -29,6 +29,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 *************************************************************************************/
 /*  END LEGAL */
 #include <Grid/Grid.h>
+#include <cinttypes>   // PRIx64: uint64_t is unsigned long on Linux, unsigned long long on macOS
 
 NAMESPACE_BEGIN(Grid);
 ///////////////////////////////////////////////////////
@@ -208,10 +209,8 @@ bool FlightRecorder::CsumLog(uint64_t hex)
         std::cerr<<"FlightRecorder::CsumLog Oops, I did it again "<< CsumLoggingCounter
 		 <<std::hex<<" "<<hex<<" "<<hexref<<std::dec<<std::endl;
 
-	fprintf(stderr,"%s:%d Oops, I did it again! Reproduce failure for csum %d %llx expect %llx\n",
-		GridHostname(),
-		GlobalSharedMemory::WorldShmRank,
-		CsumLoggingCounter,hex, hexref);
+	fprintf(stderr,"%s:%d Oops, I did it again! Reproduce failure for csum %d %" PRIx64 " expect %" PRIx64 "\n",
+		GridHostname(),GlobalSharedMemory::WorldShmRank,CsumLoggingCounter,(uint64_t)(hex),(uint64_t)(hexref));
 	BACKTRACEFP(stderr);
 	fflush(stderr);
 
@@ -302,11 +301,8 @@ void FlightRecorder::xmitLog(void *buf,uint64_t bytes)
 	fprintf(stderr,"FlightRecorder Oops step %d stage %s \n",
 		FlightRecorder::StepLoggingCounter,
 		FlightRecorder::StepName);
-	fprintf(stderr,"%s:%d Oops, send buf difference! Reproduce failure for xmit %d/%zu  %lx expect glb %lx\n",
-		GridHostname(),
-		GlobalSharedMemory::WorldShmRank,
-		XmitLoggingCounter,XmitLogVector.size(),
-		_xor, XmitLogVector[XmitLoggingCounter]); fflush(stderr);
+	fprintf(stderr,"%s:%d Oops, send buf difference! Reproduce failure for xmit %d/%zu  %" PRIx64 " expect glb %" PRIx64 "\n",
+		GridHostname(),GlobalSharedMemory::WorldShmRank,XmitLoggingCounter,XmitLogVector.size(),(uint64_t)(_xor),(uint64_t)(XmitLogVector[XmitLoggingCounter])); fflush(stderr);
 	BACKTRACEFP(stderr);
 	
 	if ( !ContinueOnFail ) GRID_ASSERT(0);
@@ -348,11 +344,8 @@ void FlightRecorder::recvLog(void *buf,uint64_t bytes,int rank)
 	fprintf(stderr,"FlightRecorder Oops step %d stage %s \n",
 		FlightRecorder::StepLoggingCounter,
 		FlightRecorder::StepName);
-	fprintf(stderr,"%s:%d Oops, recv buf difference! Reproduce failure for recv %d/%zu  %lx expect glb %lx from MPI rank %d\n",
-		GridHostname(),
-		GlobalSharedMemory::WorldShmRank,
-		RecvLoggingCounter,RecvLogVector.size(),
-		_xor, RecvLogVector[RecvLoggingCounter],rank); fflush(stderr);
+	fprintf(stderr,"%s:%d Oops, recv buf difference! Reproduce failure for recv %d/%zu  %" PRIx64 " expect glb %" PRIx64 " from MPI rank %d\n",
+		GridHostname(),GlobalSharedMemory::WorldShmRank,RecvLoggingCounter,RecvLogVector.size(),(uint64_t)(_xor),(uint64_t)(RecvLogVector[RecvLoggingCounter]),rank); fflush(stderr);
 	BACKTRACEFP(stderr);
 	
 	if ( !ContinueOnFail ) GRID_ASSERT(0);

Grid/util/Init.cc

@@ -666,7 +666,7 @@ void Grid_finalize(void)
 #endif
   std::cout<<GridLogMessage<<"******* Grid Finalize successful     ******"<<std::endl;
   std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
-  Grid_unquiesce_nodes();
+  //  Grid_unquiesce_nodes();
   Grid_is_initialised = 0;
 }
 

examples/Example_pvdagm_v2_3level_DenseCoarseMatrix.cc

@@ -82,10 +82,13 @@ std::vector<int> lat_size({48,48,48,96});
 // Solver tuning, values as in the V1 example
 RealD FineSmootherShift    = 0.1;
 int   FineSmootherOrder    = 16;
+int   FineSmootherMmax     = 1;
 RealD CoarseSmootherShift  = 0.1;
 int   CoarseSmootherNstep  = 4;
+int   CoarseSmootherMmax   = 1;
 RealD CoarseSolverTol      = 0.03;
 int   CoarseSolverOrder    = 200;
+int   CoarseSolverMmax     = 4;
 RealD OuterTol             = 1.0e-8;
 int   OuterMmax            = 8;
 int   OuterNstep           = 8;
@@ -108,6 +111,10 @@ int   OuterNstep           = 8;
 // 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)
@@ -118,13 +125,17 @@ void ParseEnvironment(void)
   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("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;
@@ -287,7 +298,7 @@ public:
   static RealD vnorm2(std::vector<Field> &x){ RealD s=0; for(auto &f:x) s+=norm2(f); return s; }
   static ComplexD vinnerProduct(std::vector<Field> &x,std::vector<Field> &y){ ComplexD s(0); for(int r=0;r<(int)x.size();r++) s+=innerProduct(x[r],y[r]); return s; }
   static void vaxpy(std::vector<Field> &z,ComplexD a,std::vector<Field> &x,std::vector<Field> &y){ for(int r=0;r<(int)z.size();r++) axpy(z[r],a,x[r],y[r]); }
-  void vOp(std::vector<Field> &in,std::vector<Field> &out){ for(int r=0;r<(int)in.size();r++) Linop.Op(in[r],out[r]); }
+  void vOp(std::vector<Field> &in,std::vector<Field> &out){ GRID_TRACE("MrhsPGCR::vOp"); for(int r=0;r<(int)in.size();r++) Linop.Op(in[r],out[r]); }
   void operator()(std::vector<Field> &src,std::vector<Field> &psi){
     RealD cp,ssq,rsq; int nrhs=src.size(); GridBase *grid=src[0].Grid();
     ssq=vnorm2(src); rsq=Tolerance*Tolerance*ssq;
@@ -324,12 +335,19 @@ public:
       vaxpy(psi,a,p[peri_k],psi); vaxpy(r,-a,q[peri_k],r); cp=vnorm2(r);
       std::cout<<GridLogMessage<<std::string(level,'\t')<<" "<<name<<" MrhsPGCR step["<<steps<<"]  resid "<<cp<<" target "<<rsq<<std::endl;
       if((k==nstep-1)||(cp<rsq)) return cp;
-      Preconditioner(r,z); vOp(z,Az); zAAz=vnorm2(Az);
+      Preconditioner(r,z);
+      vOp(z,Az);
+      zAAz=vnorm2(Az);
       q[peri_kp]=Az; p[peri_kp]=z;
       int northog=((kp)>(mmax-1))?(mmax-1):(kp);
-      for(int back=0;back<northog;back++){ int peri_back=(k-back)%mmax; GRID_ASSERT((k-back)>=0);
-        b=-real(vinnerProduct(q[peri_back],Az))/qq[peri_back];
-        vaxpy(p[peri_kp],b,p[peri_back],p[peri_kp]); vaxpy(q[peri_kp],b,q[peri_back],q[peri_kp]); }
+      {
+	GRID_TRACE("MrhsPGCR orthog");
+        for(int back=0;back<northog;back++){
+	  int peri_back=(k-back)%mmax; GRID_ASSERT((k-back)>=0);
+	  b=-real(vinnerProduct(q[peri_back],Az))/qq[peri_back];
+	  vaxpy(p[peri_kp],b,p[peri_back],p[peri_kp]); vaxpy(q[peri_kp],b,q[peri_back],q[peri_kp]);
+	}
+      }
       qq[peri_kp]=vnorm2(q[peri_kp]);
     }
     GRID_ASSERT(0); return cp;
@@ -369,14 +387,14 @@ public:
     CoarseCoarseField CCsol(_CoarseCoarseMrhs);
     _Projector.blockProject(vec1,CCsrc);
 
-    CCsol=Zero();
+    //    CCsol=Zero(); // Is this necessary?
     _CoarseCoarseSolve(CCsrc,CCsol);
 
     _Projector.blockPromote(vec1,CCsol);
     add(out,out,vec1);
 
     _CoarseOp.Op(out,vec1);  sub(vec1,in,vec1);
-    vec2=Zero();
+    //    vec2=Zero(); // Zero guess
     _CoarseSmoother(vec1,vec2);
     add(out,out,vec2);
   }
@@ -428,7 +446,10 @@ public:
       for(int r=0;r<nrhs;r++){ _FineOperator.Op(out[r],vec1[r]); sub(vec1[r],in[r],vec1[r]); }
     }
     { GRID_TRACE("MGPostSmooth");
-      for(int r=0;r<nrhs;r++){ vec2[r]=Zero(); _PostSmoother(vec1[r],vec2[r]); add(out[r],out[r],vec2[r]); }
+      for(int r=0;r<nrhs;r++){
+	//	vec2[r]=Zero();
+	_PostSmoother(vec1[r],vec2[r]); add(out[r],out[r],vec2[r]);
+      }
     }
     SetFineSloppy(0);
   }
@@ -837,7 +858,7 @@ int main (int argc, char ** argv)
 
     ShiftedLinearOperator<CoarseVector> ShiftedC(CoarseSmootherShift, LinOpC);
     PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>
-      CoarseSmootherGCR(0.01,1,ShiftedC,simpleC,CoarseSmootherNstep,CoarseSmootherNstep);
+      CoarseSmootherGCR(0.01,1,ShiftedC,simpleC,CoarseSmootherMmax,CoarseSmootherNstep);
     CoarseSmootherGCR.Level(2); CoarseSmootherGCR.Name("Csmoother"); CoarseSmootherGCR.SetZeroGuess(1);
 
     MrhsCoarseThreeLevelPrec<CoarseVector,CoarseCoarseVector>
@@ -845,11 +866,13 @@ int main (int argc, char ** argv)
                   Coarse5d, CoarseCoarse5d, CCMrhs, nr);
 
     PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>
-      L2PGCR(CoarseSolverTol, CoarseSolverOrder/16, LinOpC, L2to3Precon, 16, 16);
+      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,FineSmootherOrder,FineSmootherOrder);
+    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);
 
     MrhsTwoLevelMG<LatticeFermionD,CoarseVector,FineSmoother_t>
       ThreeLevelPrecon(PVdagM, SmootherGCR, MrhsProjector, L2PGCR, Coarse5d, CMrhs);