Stop mucking about with THREAD_MULTIPLE

benchmarks/Benchmark_allreduce.cc

@@ -48,8 +48,16 @@ Author: Peter Boyle <pboyle@bnl.gov>
 // GlobalSumVector(1) vs the two rings at n=1.
 //
 //   mpirun -n 8 ./Benchmark_allreduce --grid 16.16.16.32 --mpi 1.1.2.4
-//   env: BENCH_MIN_KB (4)  BENCH_MAX_MB (512)  BENCH_REPS (5)
+//   env: BENCH_MIN_KB (4)  BENCH_MAX_MB (512; 0 = scalar section only)  BENCH_REPS (5)
 //        BENCH_MPI_DEV_MAX_MB (4: above this MPI-dev is skipped, not attempted)
+//        BENCH_SCALAR_REPS (200)
+//        BENCH_SCALAR_GAP_US (0): a device kernel of about this many us between
+//        scalar reps.  Back-to-back sums (gap 0) measured GlobalSumP2P at 12.8 ms
+//        mean / 0.4 ms min at P=288 vs MPI 45 us (2026-08-27), yet the solver's
+//        fine-smoother steps (16 ms incl. a ~10 ms matvec) cannot be paying that
+//        per reduction: the solver spaces its sums by ms of GPU work.  The gap
+//        reproduces the solver's condition; if P2P is fast with a gap and slow
+//        without, the tight loop is the pathology, not the primitive.
 //////////////////////////////////////////////////////////////////////////////
 
 #include <Grid/Grid.h>
@@ -149,25 +157,44 @@ void Run(GridCartesian *grid, const char *tname, uint64_t nbytes_lo, uint64_t nb
   }
 }
 
-void Scalar(GridCartesian *grid, int reps)
+void Scalar(GridCartesian *grid, int reps, double gap_us)
 {
   int me=grid->ThisRank();
-  Stats<RealD> sP2P, sBare, sVec, sCart, sFlat;
+  Stats<RealD> sP2P, sBare, sVec, sCart, sFlat, sP2Pdirect;
   RealD x; deviceVector<RealD> d(1); std::vector<RealD> h(1);
+  // gap: a bandwidth-bound device kernel sized to take ~gap_us (calibrated once)
+  uint64_t gapN = 0; deviceVector<RealD> gapbuf(1);
+  if ( gap_us > 0 ) {
+    gapN = 1<<20; gapbuf.resize(gapN);
+    RealD *g=&gapbuf[0];
+    accelerator_for(i,gapN,1,{ g[i]=1.0; });
+    double t0=usecond(); accelerator_for(i,gapN,1,{ g[i]=g[i]*1.000001+1.0e-9; }); double t1=usecond();
+    double per = (t1-t0)/gapN;                              // us per element at this size
+    gapN = (uint64_t)std::max(1.0, gap_us/std::max(per,1.0e-6)); gapbuf.resize(gapN);
+    g=&gapbuf[0]; accelerator_for(i,gapN,1,{ g[i]=1.0; });
+    double t2=usecond(); accelerator_for(i,gapN,1,{ g[i]=g[i]*1.000001+1.0e-9; }); double t3=usecond();
+    if ( me==0 ) std::cout << GridLogMessage << "Scalar gap kernel: " << gapN << " elements, " << (t3-t2) << " us (requested " << gap_us << ")" << std::endl;
+  }
+  auto gap=[&](void){ if(gapN){ RealD *g=&gapbuf[0]; accelerator_for(i,gapN,1,{ g[i]=g[i]*1.000001+1.0e-9; }); } };
+  // A Barrier before each timed call so that node skew (from the gap kernel or
+  // anything else) is not charged to the reduction: the timer measures the
+  // collective from a synchronised start.
   for(int r=0;r<reps;r++){
-    x=1.0+me;   { double t0=usecond(); grid->GlobalSum(x);                 sP2P.add(usecond()-t0); }
+    gap(); x=1.0+me;   grid->Barrier(); { double t0=usecond(); grid->GlobalSum(x);                 sP2P.add(usecond()-t0); }
-    h[0]=1.0+me;{ double t0=usecond(); BareAllreduce(grid,&h[0],1);        sBare.add(usecond()-t0); }
+    gap(); x=1.0+me;   grid->Barrier(); { double t0=usecond(); grid->GlobalSumP2P(x);              sP2Pdirect.add(usecond()-t0); }
-    h[0]=1.0+me;{ double t0=usecond(); grid->GlobalSumVector(&h[0],1);     sVec.add(usecond()-t0); }
+    gap(); h[0]=1.0+me;grid->Barrier(); { double t0=usecond(); BareAllreduce(grid,&h[0],1);        sBare.add(usecond()-t0); }
-    acceleratorPut(d[0],h[0]); { double t0=usecond(); CartesianRingAllReduce(grid,&d[0],1); sCart.add(usecond()-t0); }
+    gap(); h[0]=1.0+me;grid->Barrier(); { double t0=usecond(); grid->GlobalSumVector(&h[0],1);     sVec.add(usecond()-t0); }
-    acceleratorPut(d[0],h[0]); { double t0=usecond(); RingAllReduce(grid,&d[0],1);          sFlat.add(usecond()-t0); }
+    gap(); acceleratorPut(d[0],h[0]); grid->Barrier(); { double t0=usecond(); CartesianRingAllReduce(grid,&d[0],1); sCart.add(usecond()-t0); }
+    gap(); acceleratorPut(d[0],h[0]); grid->Barrier(); { double t0=usecond(); RingAllReduce(grid,&d[0],1);          sFlat.add(usecond()-t0); }
   }
   auto mx=[&](double v){ RealD y=v; grid->GlobalMax(y); return (double)y; };
   auto line=[&](const char *nm, Stats<RealD> &st){
     double tmin=mx(st.tmin), tmean=mx(st.tsum/st.n);
     if ( me==0 ) std::cout << GridLogMessage << "  " << std::setw(34) << std::left << nm << std::right
                            << " min " << std::setw(8) << tmin << " us   mean " << std::setw(8) << tmean << " us" << std::endl; };
-  if ( me==0 ) std::cout << GridLogMessage << "==== SCALAR latency (RealD, " << reps << " reps, slowest rank)  P=" << grid->ProcessorCount() << std::endl;
+  if ( me==0 ) std::cout << GridLogMessage << "==== SCALAR latency (RealD, " << reps << " reps, slowest rank, Barrier before each timed call, gap " << gap_us << " us)  P=" << grid->ProcessorCount() << std::endl;
   line("GlobalSum(RealD) = GlobalSumP2P",   sP2P);
+  line("GlobalSum(RealD) = GlobalSumP2Pdirect",   sP2Pdirect);
   line("bare MPI_Allreduce(1 double)",       sBare);
   line("GlobalSumVector(double*,1) [MPI]",   sVec);
   line("CartesianRingAllReduce n=1 (device)",sCart);
@@ -178,7 +205,8 @@ int main(int argc, char **argv)
 {
   Grid_init(&argc, &argv);
   GridCartesian *grid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd, vComplexD::Nsimd()), GridDefaultMpi());
-  Scalar(grid, getenv("BENCH_SCALAR_REPS") ? atoi(getenv("BENCH_SCALAR_REPS")) : 200);
+  Scalar(grid, getenv("BENCH_SCALAR_REPS") ? atoi(getenv("BENCH_SCALAR_REPS")) : 200,
+               getenv("BENCH_SCALAR_GAP_US") ? atof(getenv("BENCH_SCALAR_GAP_US")) : 0.0);
   uint64_t lo   = (getenv("BENCH_MIN_KB")         ? atol(getenv("BENCH_MIN_KB"))         : 4) * 1024ull;
   uint64_t hi   = (getenv("BENCH_MAX_MB")         ? atol(getenv("BENCH_MAX_MB"))         : 512) * 1024ull*1024ull;
   int      reps = getenv("BENCH_REPS")            ? atoi(getenv("BENCH_REPS"))            : 5;

examples/Example_pvdagm_v2_3level_DenseCoarseMatrix.cc

@@ -556,17 +556,6 @@ 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();
 

systems/Frontier/allreduce_bench.job

@@ -4,7 +4,7 @@
 #SBATCH --ntasks-per-node=8
 #SBATCH --cpus-per-task=7
 #SBATCH --gpus-per-node=8
-#SBATCH --time=0:30:00
+#SBATCH --time=0:20:00
 #SBATCH --account=phy157_dwf
 #SBATCH --gpu-bind=none
 #SBATCH --exclusive
@@ -66,9 +66,29 @@ run_cell () {
 }
 
 #          name        nodes ntask  geom      vol
-run_cell   N1_8gcd     1     8      1.1.2.4   16.16.32.64
+#run_cell   N1_8gcd     1     8      1.1.2.4   16.16.32.64      # done 2026-08-27
-run_cell   N4_32gcd    4     32     2.2.2.4   32.32.32.64
+#run_cell   N4_32gcd    4     32     2.2.2.4   32.32.32.64      # done: MPI <1 MB, rings above, 2-2.5x at 512 MB
-run_cell   N36_288gcd  36    288    3.6.4.4   48.48.48.96
+#run_cell   N36_288gcd  36    288    3.6.4.4   48.48.48.96      # done: MPI <2 MB, cart 2-48 MB, flat above (3x at 512 MB);
+                                                               #       SCALAR: GlobalSumP2P 0.4 ms min / 12.8 ms MEAN vs MPI 45 us (!)
+
+# SCALAR-ONLY sweep at 288 ranks (BENCH_MAX_MB=0 skips the vector tables; each cell ~30 s):
+# S0 back-to-back as measured; S1 with ~2 ms of GPU work between sums (the solver's condition:
+# fine-smoother steps are 16 ms apart incl. a ~10 ms matvec, so the solver cannot be paying
+# 12 ms per reduction); S2-S3 one environment knob each.  Readout: the five "min/mean us" lines.
+export BENCH_MAX_MB=0
+scalar_cell () { name=$1; shift; run_cell_env "$name" 36 288 3.6.4.4 48.48.48.96 "$@"; }
+run_cell_env () {
+  name=$1; nodes=$2; ntask=$3; geom=$4; vol=$5; shift 5
+  echo "----- $name : $* -----"
+  fname=log.allreduce.$name
+  env "$@" srun -N$nodes -n$ntask --kill-on-bad-exit=1 ./select_gpu $BIN --mpi $geom --grid $vol $OPTS1 > $fname 2>&1
+  echo "  exit $?"; sleep 20
+  grep -h "==== SCALAR\|gap kernel\| min .* us" $fname | sed 's/^Grid : Message : [0-9.]* s : //' | cut -c1-140
+}
+scalar_cell  S0_tight        BENCH_SCALAR_GAP_US=0
+scalar_cell  S1_gap2ms       BENCH_SCALAR_GAP_US=2000
+scalar_cell  S2_gap2ms_omp1  BENCH_SCALAR_GAP_US=2000 OMP_NUM_THREADS=1
+scalar_cell  S3_gap2ms_nogpu BENCH_SCALAR_GAP_US=2000 MPICH_GPU_SUPPORT_ENABLED=0
 
 echo "========================================================="
 echo "done: see log.allreduce.N1_8gcd, N4_32gcd, N36_288gcd"

systems/Frontier/schur2d_env.job

@@ -76,7 +76,7 @@ run_cell () {
 
 #          name              environment
 run_cell   E1_omp1           OMP_NUM_THREADS=1
-run_cell   E2_omp1_multiple  OMP_NUM_THREADS=1 GRID_MPI_THREAD_MULTIPLE=1 MPICH_MAX_THREAD_SAFETY=multiple     # = the SLATE harness environment
+#run_cell  E2_omp1_multiple  (MPI_THREAD_MULTIPLE init: 21.2 s, NIL; the GRID_MPI_THREAD_MULTIPLE hack was removed 2026-08-28)
 run_cell   E3_omp7           OMP_NUM_THREADS=7                                                                 # the example's thread count
 run_cell   E4_omp7_ballast   OMP_NUM_THREADS=7 S2D_BALLAST_GB=10.6                                            # + the example's device residency
 run_cell   E5_omp1_ballast   OMP_NUM_THREADS=1 S2D_BALLAST_GB=10.6                                            # ballast alone

systems/Frontier/smoother_modes.job

@@ -102,10 +102,8 @@ 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);
+# Inverse ring-rate hypotheses (2026-08-27): OMP threads and MPI thread level both NIL
-# (2) if (1) fails, uncomment the two lines below (harness ran 62 s with these).
+# (schur2d_env.job); the cause was partial ring participation -> SCHUR2D_LEAF_SPAN.
-#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

tests/debug/Test_schur2d_scale.cc

@@ -76,20 +76,13 @@ static ComplexD Fill(int64_t i, int64_t j, int64_t N)
 
 int main(int argc, char **argv)
 {
-  // Environment walk (2026-08-27): the SAME inverse runs 20.8 s in the SLATE
+  // Environment walk (2026-08-27, systems/Frontier/schur2d_env.job): knobs to
-  // harness and 26.5 s in the example, with the local GPU work (GEMM, leaf
+  // reproduce the example's environment here, one at a time.  Result: thread
-  // inverse) 25-40% slower in the example.  Knobs to reproduce the example's
+  // level, OMP_NUM_THREADS, device residency and sustained load all NIL; only
-  // environment here, one at a time (systems/Frontier/schur2d_env.job):
+  // "first job step on fresh nodes" (+3 s) is real.
-  //   GRID_MPI_THREAD_MULTIPLE=1  MPI_Init_thread(MULTIPLE) before Grid_init (SLATE harness does this)
+  //   S2D_BALLAST_GB=x     x GB of Lattice fields made device-resident before the invert
-  //   S2D_BALLAST_GB=x            x GB of Lattice fields made device-resident before the invert
+  //   S2D_PREHEAT_S=x      x seconds of back-to-back zgemm before the invert
-  //                               (the example carries ~10.6 GB of fine-grid state in the MemoryManager)
+  //   OMP_NUM_THREADS      set in the job, read by nothing here but the runtime
-  //   OMP_NUM_THREADS             set in the job, read by nothing here but the runtime
-  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);
 
   GridCartesian *grid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(),