More tracing

Grid/algorithms/deflation/MultiRHSBlockProject.h

@@ -760,6 +760,7 @@ public:
   template<class cobj>
   void blockProject(std::vector<Field> &fine,std::vector< Lattice<cobj> > & coarse)
   {
+    GRID_TRACE("BlockProject");
     int nrhs=fine.size();
     int _nbasis = sizeof(typename cobj::scalar_object)/sizeof(scalar);
     //    std::cout << "blockProject nbasis " <<nbasis<<" " << _nbasis<<std::endl;
@@ -819,6 +820,7 @@ public:
   template<class cobj>
   void blockPromote(std::vector<Field> &fine,std::vector<Lattice<cobj> > & coarse)
   {
+    GRID_TRACE("BlockPromote");
     int nrhs=fine.size();
     int _nbasis = sizeof(typename cobj::scalar_object)/sizeof(scalar);
     GRID_ASSERT(nbasis==_nbasis);
@@ -877,6 +879,7 @@ public:
   template<class cobj>
   void blockProject(Field &fine_mrhs,Lattice<cobj> &coarse_mrhs)
   {
+    GRID_TRACE("BlockProjectMrhs");
     int nrhs = fine_mrhs.Grid()->_rdimensions[0];
     int _nbasis = sizeof(typename cobj::scalar_object)/sizeof(scalar);
     GRID_ASSERT(nbasis==_nbasis);
@@ -893,6 +896,7 @@ public:
   template<class cobj>
   void blockPromote(Field &fine_mrhs,Lattice<cobj> &coarse_mrhs)
   {
+    GRID_TRACE("BlockPromoteMrhs");
     int nrhs = fine_mrhs.Grid()->_rdimensions[0];
     int _nbasis = sizeof(typename cobj::scalar_object)/sizeof(scalar);
     GRID_ASSERT(nbasis==_nbasis);
@@ -913,6 +917,7 @@ public:
   template<class cobj>
   void blockProject(std::vector<Field> &fine,Lattice<cobj> &coarse_mrhs)
   {
+    GRID_TRACE("BlockProjectMixed");
     int nrhs = fine.size();
     int _nbasis = sizeof(typename cobj::scalar_object)/sizeof(scalar);
     GRID_ASSERT(nbasis==_nbasis);
@@ -929,6 +934,7 @@ public:
   template<class cobj>
   void blockProject(Field &fine_mrhs,std::vector< Lattice<cobj> > &coarse)
   {
+    GRID_TRACE("BlockProjectMixed");
     int nrhs = fine_mrhs.Grid()->_rdimensions[0];
     int _nbasis = sizeof(typename cobj::scalar_object)/sizeof(scalar);
     GRID_ASSERT(nbasis==_nbasis);
@@ -945,6 +951,7 @@ public:
   template<class cobj>
   void blockPromote(std::vector<Field> &fine,Lattice<cobj> &coarse_mrhs)
   {
+    GRID_TRACE("BlockPromoteMixed");
     int nrhs = fine.size();
     int _nbasis = sizeof(typename cobj::scalar_object)/sizeof(scalar);
     GRID_ASSERT(nbasis==_nbasis);
@@ -961,6 +968,7 @@ public:
   template<class cobj>
   void blockPromote(Field &fine_mrhs,std::vector< Lattice<cobj> > &coarse)
   {
+    GRID_TRACE("BlockPromoteMixed");
     int nrhs = fine_mrhs.Grid()->_rdimensions[0];
     int _nbasis = sizeof(typename cobj::scalar_object)/sizeof(scalar);
     GRID_ASSERT(nbasis==_nbasis);
@@ -1002,6 +1010,7 @@ public:
   // C_br = V^dag F
   void ProjectBLAS(int nrhs)
   {
+    GRID_TRACE("ProjectBLAS");
     deviceVector<scalar *> Vd(coarse_vol);
     deviceVector<scalar *> Fd(coarse_vol);
     deviceVector<scalar *> Cd(coarse_vol);
@@ -1022,6 +1031,7 @@ public:
   // F_xr = Vxb Cbr
   void PromoteBLAS(int nrhs)
   {
+    GRID_TRACE("PromoteBLAS");
     deviceVector<scalar *> Vd(coarse_vol);
     deviceVector<scalar *> Fd(coarse_vol);
     deviceVector<scalar *> Cd(coarse_vol);

Grid/algorithms/multigrid/DenseCoarseMatrix.h

@@ -948,6 +948,7 @@ public:
   ////////////////////////////////////////////////////////////////////
   void SlabApplyPacked(int nr, double *tprof)
   {
+    GRID_TRACE("DenseSlabApply");
     GRID_ASSERT(nr <= MRHS_MAX);
     uint64_t nX = (uint64_t)N * nr;
     uint64_t nY = (uint64_t)nrows * nr;
@@ -955,14 +956,22 @@ public:
     double t1 = usecond();
     double t2, t3;
     if (devSum) {
-      acceleratorCopyToDevice(&hX[0],&dX[0],nX*sizeof(ComplexF));
+      { GRID_TRACE("DenseH2D");
+        acceleratorCopyToDevice(&hX[0],&dX[0],nX*sizeof(ComplexF));
+      }
       t2 = usecond();
-      grid->GlobalSumVector((ComplexF *)&dX[0], (int)nX);
+      { GRID_TRACE("DenseAllreduce");
+        grid->GlobalSumVector((ComplexF *)&dX[0], (int)nX);
+      }
       t3 = usecond();
     } else {
-      grid->GlobalSumVector(&hX[0], (int)nX);
+      { GRID_TRACE("DenseAllreduce");
+        grid->GlobalSumVector(&hX[0], (int)nX);
+      }
       t2 = usecond();
-      acceleratorCopyToDevice(&hX[0],&dX[0],nX*sizeof(ComplexF));
+      { GRID_TRACE("DenseH2D");
+        acceleratorCopyToDevice(&hX[0],&dX[0],nX*sizeof(ComplexF));
+      }
       t3 = usecond();
     }
     // Y = op(slab,T) . X : row-major slab (nrows x N) == col-major A^T
@@ -970,13 +979,13 @@ public:
     // Split-K: NK chunk-GEMMs by pointer offset (AOT lists), then reduce.
     ComplexF one (1.0,0.0);
     ComplexF zero(0.0,0.0);
-    BLAS.gemmBatched(GridBLAS_OP_T, GridBLAS_OP_N,
+    { GRID_TRACE("DenseSplitKGEMM");
-                     (int)nrows, nr, (int)Kc,
+      BLAS.gemmBatched(GridBLAS_OP_T, GridBLAS_OP_N,
-                     one,  aptrs, (int)N,
+                       (int)nrows, nr, (int)Kc,
-                           xptrs, (int)N,
+                       one,  aptrs, (int)N,
-                     zero, cptrs, (int)nrows);
+                             xptrs, (int)N,
-    BLAS.synchronise();
+                       zero, cptrs, (int)nrows);
-    {
+      BLAS.synchronise();
       ComplexF *pp = &dPartial[0];
       ComplexF *py = &dY[0];
       uint64_t stride = (uint64_t)nrows*MRHS_MAX;
@@ -988,7 +997,9 @@ public:
       });
     }
     double t4 = usecond();
-    acceleratorCopyFromDevice(&dY[0],&hY[0],nY*sizeof(ComplexF));
+    { GRID_TRACE("DenseD2H");
+      acceleratorCopyFromDevice(&dY[0],&hY[0],nY*sizeof(ComplexF));
+    }
     double t5 = usecond();
     if (tprof) {
       tprof[0] = devSum ? (t3-t2) : (t2-t1);   // allreduce
@@ -1003,20 +1014,25 @@ public:
   ////////////////////////////////////////////////////////////////////
   virtual void operator()(const Field &src, Field &psi)
   {
+    GRID_TRACE("DenseApply1");
     uint64_t nX = (uint64_t)N;
-    thread_for(i, nX, { hX[i]=ComplexF(0.0,0.0); });
+    { GRID_TRACE("DensePack");
-    for(int ss=0; ss<lsites; ss++){
+      thread_for(i, nX, { hX[i]=ComplexF(0.0,0.0); });
-      sobj s;
+      for(int ss=0; ss<lsites; ss++){
-      peekLocalSite(s, src, myLcoor[ss]);
+        sobj s;
-      for(int b=0; b<nbasis; b++)
+        peekLocalSite(s, src, myLcoor[ss]);
-        hX[ myGsite[ss]*nbasis + b ] = ComplexF(((ComplexD *)&s)[b]);
+        for(int b=0; b<nbasis; b++)
+          hX[ myGsite[ss]*nbasis + b ] = ComplexF(((ComplexD *)&s)[b]);
+      }
     }
     SlabApplyPacked(1, nullptr);
-    for(int ss=0; ss<lsites; ss++){
+    { GRID_TRACE("DenseUnpack");
-      sobj s;
+      for(int ss=0; ss<lsites; ss++){
-      for(int b=0; b<nbasis; b++)
+        sobj s;
-        ((ComplexD *)&s)[b] = ComplexD(hY[ss*nbasis + b]);
+        for(int b=0; b<nbasis; b++)
-      pokeLocalSite(s, psi, myLcoor[ss]);
+          ((ComplexD *)&s)[b] = ComplexD(hY[ss*nbasis + b]);
+        pokeLocalSite(s, psi, myLcoor[ss]);
+      }
     }
   }
 
@@ -1044,22 +1060,26 @@ public:
     GRID_ASSERT(nr <= MRHS_MAX);
     double t0 = usecond();
     uint64_t nX = (uint64_t)N*nr;
-    thread_for(i, nX, { hX[i]=ComplexF(0.0,0.0); });
+    { GRID_TRACE("DensePack");
-    for(int rr=0; rr<nr; rr++){
+      thread_for(i, nX, { hX[i]=ComplexF(0.0,0.0); });
-      for(int ss=0; ss<lsites; ss++){
+      for(int rr=0; rr<nr; rr++){
-        sobj s;
+        for(int ss=0; ss<lsites; ss++){
-        peekLocalSite(s, src[rr], myLcoor[ss]);
+          sobj s;
-        for(int b=0; b<nbasis; b++)
+          peekLocalSite(s, src[rr], myLcoor[ss]);
-          hX[ (uint64_t)rr*N + myGsite[ss]*nbasis + b ] = ComplexF(((ComplexD *)&s)[b]);
+          for(int b=0; b<nbasis; b++)
+            hX[ (uint64_t)rr*N + myGsite[ss]*nbasis + b ] = ComplexF(((ComplexD *)&s)[b]);
+        }
       }
     }
     SlabApplyPacked(nr, nullptr);
-    for(int rr=0; rr<nr; rr++){
+    { GRID_TRACE("DenseUnpack");
-      for(int ss=0; ss<lsites; ss++){
+      for(int rr=0; rr<nr; rr++){
-        sobj s;
+        for(int ss=0; ss<lsites; ss++){
-        for(int b=0; b<nbasis; b++)
+          sobj s;
-          ((ComplexD *)&s)[b] = ComplexD(hY[(uint64_t)rr*nrows + (ss*nbasis+b)]);
+          for(int b=0; b<nbasis; b++)
-        pokeLocalSite(s, psi[rr], myLcoor[ss]);
+            ((ComplexD *)&s)[b] = ComplexD(hY[(uint64_t)rr*nrows + (ss*nbasis+b)]);
+          pokeLocalSite(s, psi[rr], myLcoor[ss]);
+        }
       }
     }
     double t1 = usecond();
@@ -1090,7 +1110,7 @@ public:
     Field &in = const_cast<Field &>(in6);
     uint64_t nX = (uint64_t)N * nr;
     thread_for(i, nX, { hX[i]=ComplexF(0.0,0.0); });
-    {
+    { GRID_TRACE("DensePack");
       autoView(iv, in, CpuRead);
       Coordinate c6(nd+1);
       for(int ss=0; ss<lsites; ss++){
@@ -1108,7 +1128,7 @@ public:
     double tprof[4];
     SlabApplyPacked(nr, tprof);
     double t5 = usecond();
-    {
+    { GRID_TRACE("DenseUnpack");
       autoView(ov, out6, CpuWrite);
       Coordinate c6(nd+1);
       for(int ss=0; ss<lsites; ss++){

Grid/algorithms/multigrid/GeneralCoarsenedMatrixMultiRHSV2.h

@@ -785,10 +785,14 @@ public:
       GRID_ASSERT(in.Grid() == _CoarseGridMulti);
     }
 
+    GRID_TRACE("CoarseV2Mult");
     t_tot=-usecond();
     CoarseVector tin=in;
     t_exch=-usecond();
-    CoarseVector pin = CellMulti->ExchangePeriodic(tin); //padded input
+    // lambda scope so the roctx range covers exactly the exchange; the
+    // PaddedCellFwd/BwdMPI markers inside it then nest properly.
+    CoarseVector pin = [&](){ GRID_TRACE("CoarseV2Exchange");
+                              return CellMulti->ExchangePeriodic(tin); }(); //padded input
     t_exch+=usecond();
 
     CoarseVector pout(pin.Grid());
@@ -812,12 +816,15 @@ public:
     
 
     t_GtoB=-usecond();
-    GridtoBLAS(pin,BLAS_B);
+    { GRID_TRACE("CoarseV2GridToBLAS");
+      GridtoBLAS(pin,BLAS_B);
+    }
     t_GtoB+=usecond();
 
     GridBLAS BLAS;
 
     t_mult=-usecond();
+    { GRID_TRACE("CoarseV2StencilGEMM");
     for(int p=0;p<geom.npoint;p++){
       RealD c = 1.0;
       if (p==0) c = 0.0;
@@ -831,10 +838,13 @@ public:
 		       BLAS_CP);
     }
     BLAS.synchronise();
+    }
     t_mult+=usecond();
 
     t_BtoG=-usecond();
-    BLAStoGrid(out,BLAS_C);
+    { GRID_TRACE("CoarseV2BLASToGrid");
+      BLAStoGrid(out,BLAS_C);
+    }
     t_BtoG+=usecond();
     t_tot+=usecond();
     /*

examples/Example_pvdagm_v2_3level_DenseCoarseMatrix.cc

@@ -90,13 +90,25 @@ RealD OuterTol             = 1.0e-8;
 int   OuterMmax            = 8;
 int   OuterNstep           = 8;
 
-// Halo exchange in reduced precision on the fine stencils.  Stencil::SloppyComms
+// "It's legal to get the same answer faster, not to get a less correct
-// is a plain runtime setter (Stencil.h:303) applied to the full/even/odd stencils
+// answer."  (PB, 2026-08-24)
-// of an operator that is already built, so this costs no extra storage and no
+//
-// second operator.  Benchmark_dwf at this decomposition measures ~1.6x on the
+// Halo-precision POLICY: reduced-precision (fp32 wire)
-// fine Dhop, consistent with the uncompressed halo being ~1.6x the interior
+// halos belong in the PRECONDITIONER -- the smoother, the V-cycle's own
-// compute time.  Set FineSloppyComms=0 to recover the exact-comms behaviour.
+// 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;
+std::function<void(int)> SetFineSloppy = [](int){};
 
 void ParseEnvironment(void)
 {
@@ -388,23 +400,37 @@ public:
     : _FineOperator(FineOp),_PostSmoother(Post),_Projector(Projector),_CoarseSolve(CoarseSolve),
       _CoarseGrid(CoarseGrid),_CoarseGridMrhs(CoarseGridMrhs){}
   virtual void operator()(std::vector<FineField> &in, std::vector<FineField> &out){
+    // The whole V-cycle is preconditioner: its fine residuals and the
+    // smoother run with sloppy halos; the caller (the outer Krylov) gets
+    // the exact operator back on exit.
+    GRID_TRACE("MGVcycle");
+    SetFineSloppy(FineSloppyComms);
     int nrhs=in.size(); GridBase *fgrid=in[0].Grid();
     std::vector<FineField> vec1(nrhs,fgrid),vec2(nrhs,fgrid);
     for(int r=0;r<nrhs;r++) out[r]=in[r];
-    for(int r=0;r<nrhs;r++){ _FineOperator.Op(out[r],vec1[r]); sub(vec1[r],in[r],vec1[r]); }
+    { GRID_TRACE("MGFineResidual");
-
+      for(int r=0;r<nrhs;r++){ _FineOperator.Op(out[r],vec1[r]); sub(vec1[r],in[r],vec1[r]); }
+    }
     // fine vector -> D+1 coarse, via the mixed blockProject
     CoarseVector CsrcMrhs(_CoarseGridMrhs), CsolMrhs(_CoarseGridMrhs);
-    _Projector.blockProject(vec1,CsrcMrhs);
+    { GRID_TRACE("MGProject");
-
+      _Projector.blockProject(vec1,CsrcMrhs);
+    }
     CsolMrhs=Zero();
-    _CoarseSolve(CsrcMrhs,CsolMrhs);
+    { GRID_TRACE("MGCoarseSolve");
-
+      _CoarseSolve(CsrcMrhs,CsolMrhs);
-    _Projector.blockPromote(vec1,CsolMrhs);
+    }
-    for(int r=0;r<nrhs;r++) add(out[r],out[r],vec1[r]);
+    { GRID_TRACE("MGPromote");
-
+      _Projector.blockPromote(vec1,CsolMrhs);
-    for(int r=0;r<nrhs;r++){ _FineOperator.Op(out[r],vec1[r]); sub(vec1[r],in[r],vec1[r]); }
+      for(int r=0;r<nrhs;r++) add(out[r],out[r],vec1[r]);
-    for(int r=0;r<nrhs;r++){ vec2[r]=Zero(); _PostSmoother(vec1[r],vec2[r]); add(out[r],out[r],vec2[r]); }
+    }
+    { GRID_TRACE("MGFineResidual2");
+      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]); }
+    }
+    SetFineSloppy(0);
   }
 };
 
@@ -484,15 +510,18 @@ int main (int argc, char ** argv)
   MobiusFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c);
   MobiusFermionD Dpv (Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,1.0, M5,b,c);
 
-  // Reduced-precision halo on both fine operators.  PVdagM and ShiftedPVdagM are
+  // PVdagM and ShiftedPVdagM are thin wrappers over these same two objects,
-  // thin wrappers over these same objects, so this covers coarsening, the fine
+  // so this one callback controls every fine halo in the program.  Default
-  // smoother, the V-cycle and the outer Krylov alike.  The final residual check
+  // EXACT; the preconditioner and the coarsening turn sloppiness on for
-  // below turns it off again: a verification computed with a sloppy operator
+  // their own scope only (policy note at FineSloppyComms).
-  // would certify the wrong matrix.
+  SetFineSloppy = [&Ddwf,&Dpv](int sloppy){
-  Ddwf.SloppyComms(FineSloppyComms);
+    Ddwf.SloppyComms(sloppy);
-  Dpv .SloppyComms(FineSloppyComms);
+    Dpv .SloppyComms(sloppy);
-  std::cout << GridLogMessage << "Fine stencils: SloppyComms = " << FineSloppyComms
+  };
-            << (FineSloppyComms ? "  (reduced-precision halo)" : "  (exact halo)") << std::endl;
+  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;
@@ -546,6 +575,7 @@ int main (int argc, char ** argv)
   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
@@ -556,6 +586,7 @@ int main (int argc, char ** argv)
     // PVdagM is single RHS: apply it directly, batch on the coarse side.
     CoarseOpPV.CoarsenOperator(PVdagM,AggregatesGCR.subspace,Coarse5d,batch);
   }
+  SetFineSloppy(0);
 
   // 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.
@@ -602,7 +633,9 @@ int main (int argc, char ** argv)
 
     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);
@@ -834,10 +867,9 @@ int main (int argc, char ** argv)
     std::cout << GridLogMessage << "V2 3-level solve Nrhs "<<nr<<" total " << w.Elapsed()
               << "  (per RHS: " << w.useconds()/1.0e6/nr << " s)" << std::endl;
 
-    // Verify against the EXACT operator: the solve may have used a reduced
+    // The outer operator is exact by policy; assert the state rather than
-    // precision halo, but the residual we report must not.
+    // trust it -- a preconditioner that failed to restore would surface here.
-    Ddwf.SloppyComms(0);
+    SetFineSloppy(0);
-    Dpv .SloppyComms(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]));
@@ -846,8 +878,7 @@ int main (int argc, char ** argv)
       std::cout << GridLogMessage << "FINAL Nrhs "<<nr<<": worst-case residual = " << worst
                 << "   (exact-halo verification)" << std::endl;
     }
-    Ddwf.SloppyComms(FineSloppyComms);
+
-    Dpv .SloppyComms(FineSloppyComms);
 
     // The operators borrow these grids and build a PaddedCell on them, so
     // they must let go before the grids are destroyed.