Test script for Frontier on I/O

Switch on the faster IO be default
IO updates for AllToAllV aggregation of contiguous chunks.
2026-11-02 18:28:04 +00:00 · 2026-08-11 13:37:50 -04:00 · 2026-08-11 13:36:50 -04:00 · 2026-08-11 13:33:04 -04:00 · 2026-06-26 11:50:44 -04:00 · 2026-06-17 20:47:15 +02:00
34 changed files with 2004 additions and 1182 deletions
@@ -30,9 +30,6 @@ Key configure options:
 | `--enable-Nc=` | `3` (default), `2`, `4`, `5` |
 | `--with-gmp=`, `--with-mpfr=`, `--with-fftw=`, `--with-lime=` | paths to libs |
 | `--enable-hdf5`, `--enable-mkl`, `--enable-lapack` | optional features |
-| `--enable-debug=yes` | adds `-g`, removes `-O3` |
-
-Use `make V=1` for verbose compiler output (shows full flags; required for bug reports).

 Platform recipes from `README.md`:
 - **KNL**: `--enable-simd=KNL --enable-comms=mpi3-auto --enable-mkl`
@@ -42,17 +39,39 @@ Platform recipes from `README.md`:

 Required external libs: GMP, MPFR, OpenSSL, zlib.

-## Running Tests
+### Use `systems/` for real machines
+
+`systems/<machine>/` holds the known-good build for each production platform (`Frontier`, `Aurora`, `Perlmutter`, `Summit`, `Tursa`, `Lumi`, `Booster`, `Crusher`, `SDCC-*`, `mac-arm`, …). Each contains a `config-command` (the exact `../../configure` invocation) and a `sourceme.sh` (module loads and env). **Prefer copying/adapting these over hand-rolling configure flags** — they encode compiler workarounds, `LDFLAGS`, and shared-memory settings that are easy to get wrong. `systems/WorkArounds.txt` records known vendor bugs.
+
+Note the GPU builds use `--enable-simd=GPU --enable-gen-simd-width=64`, so `Nsimd` is *not* 1 on device (it is `64/sizeof(scalar)`).
+
+### Regenerating `Make.inc` — required after adding or deleting source files
+
+`Make.inc` files are generated, not tracked in git (`.gitignore`d). `scripts/filelist` walks `Grid/`, `tests/*`, `benchmarks/`, `examples/`, and `HMC/` and writes the file lists and per-test `bin_PROGRAMS` rules. Every new `.cc`/`.h` in `Grid/`, and every new `Test_*.cc` / `Benchmark_*.cc` / `Example_*.cc`, is invisible to the build until you run:
+
+```bash
+./scripts/filelist    # from the source root, then re-run configure/make
+```
+
+`bootstrap.sh` runs it for you on the first setup.
+
+## Running Tests and Benchmarks

 ```bash
 # From build directory
-make check                          # root-level tests (Test_simd, Test_cshift, etc.)
-make -C tests/<subdir> tests        # build tests in a subdirectory
+make check                          # runs only Test_simd, Test_cshift, Test_stencil, Test_dwf_mixedcg_prec
+make -C tests/<subdir> tests        # build (not run) the tests in a subdirectory
+make bench                          # build benchmarks
 ./tests/core/Test_simd              # run a single test binary directly
+mpirun -n 4 ./tests/core/Test_cshift --grid 16.16.16.16 --mpi 1.1.1.4
 ```

+`make check` is a thin smoke test — building a subdirectory with `make -C tests/<subdir> tests` and running the relevant binaries directly is the normal development loop. Test binaries take Grid's standard command-line arguments (`--grid`, `--mpi`, `--accelerator-threads`, `--threads`, `--debug-signals`, `--log`); see `Grid/util/Init.cc`.
+
 Test subdirectories and their focus: `core` (SIMD, stencil, comms), `solver` (CG, GMRES, eigensolvers), `hmc` (MD integrators), `forces` (fermion forces), `lanczos`, `IO`, `smearing`, `sp2n`, `debug`.

+Tests and benchmarks that need optional fermion representations are guarded by `disable_tests_without_instantiations.h` / `disable_benchmarks_without_instantiations.h`, so a `--disable-fermion-reps --disable-gparity` build silently compiles them to no-ops.
+
 ## Architecture

 ### Layer stack (bottom to top)
@@ -78,9 +97,35 @@ Test subdirectories and their focus: `core` (SIMD, stencil, comms), `solver` (CG
 - `smearing/` — APE, Stout, HEX, gradient flow
 - `observables/` — Polyakov loop, plaquette, topological charge

-### GPU acceleration
+### GPU acceleration and the view/memory-manager discipline

-GPU support is injected via macros (`accelerator_for`, `accelerator_for2dNB`). The `Grid/simd/` SIMD types map to scalar on GPU device code; host code paths remain vectorised. Unified virtual memory is on by default (`--enable-unified=yes`); device-aware MPI (`--enable-accelerator-aware-mpi`) avoids device→host copies on transfers.
+GPU support is injected via macros in `Grid/threads/Accelerator.h` — `accelerator_for(i, n, nsimd, {...})`, `accelerator_forNB` (non-blocking, must be followed by `accelerator_barrier()`), `accelerator_for2dNB`, and `accelerator_inline`. On a CPU build these degrade to `thread_for` (OpenMP). Unified virtual memory is on by default (`--enable-unified=yes`); device-aware MPI (`--enable-accelerator-aware-mpi`) avoids device→host copies on transfers.
+
+Lattice data is **not** directly addressable inside a kernel. You must open a view with the correct access mode so `Grid/allocator/MemoryManager.h` can move/mark the data:
+
+```cpp
+autoView(out_v, out, AcceleratorWriteDiscard);   // RAII; closes at end of scope
+autoView(in_v,  in,  AcceleratorRead);
+accelerator_for(ss, grid->oSites(), Nsimd, {
+    coalescedWrite(out_v[ss], coalescedRead(in_v[ss]));
+});
+```
+
+Modes are `AcceleratorRead/Write/WriteDiscard` and `CpuRead/Write/WriteDiscard`. Getting the mode wrong (e.g. `AcceleratorRead` on a field you write) produces stale-data bugs that only appear on GPU builds. Inside kernels use `coalescedRead`/`coalescedWrite` rather than raw `operator[]` — they map the SIMD lane onto `threadIdx.x` so accesses stay coalesced.
+
+### Repo-local debugging skills (`skills/`)
+
+`skills/` contains hard-won, Grid-specific playbooks written as invocable skill files. Consult them before debugging in these areas rather than reasoning from first principles:
+
+| File | Covers |
+|---|---|
+| `gpu-memory-performance.md` | `acceleratorThreads()`, LambdaApply thread mapping, `coalescedRead` idiom, fused vs staged HBM access |
+| `gpu-runtime-correctness.md` | GPU runtime returning early from sync, silent wrong answers |
+| `communication-overlap.md` | 7-phase halo pipeline, per-packet events, host-staging vs GPU-direct RDMA |
+| `mpi-heterogeneous.md` | `MPI_Sendrecv` device-buffer aliasing, deterministic reductions |
+| `compiler-validation.md` | Isolating GPU compiler codegen bugs, minimal reproducers |
+| `correctness-verification.md` | Double-run fingerprinting, per-packet checksums, flight recorder |
+| `hang-diagnosis.md` | Diagnosing MPI/accelerator hangs |

 ### Memory and I/O

@@ -88,28 +133,24 @@ GPU support is injected via macros (`accelerator_for`, `accelerator_for2dNB`). T
 - `Grid/parallelIO/` — distributed parallel reader/writer for ILDG (via LIME), SciDAC, and native binary formats
 - `Grid/serialisation/` — text, binary, HDF5, XML/JSON serialisation of arbitrary Grid objects

-### HMC applications
+### Executables

-`HMC/` contains production-ready HMC driver programmes (e.g. `Mobius2p1f.cc`, `DWF_plus_DSDR_nf2plus1_Shamir_Gparity.cc`). These are built separately from the library tests.
+- `HMC/` — production HMC driver programmes (e.g. `Mobius2p1f.cc`, `DWF_plus_DSDR_nf2plus1_Shamir_Gparity.cc`)
+- `benchmarks/` — `Benchmark_dwf`, `Benchmark_ITT`, `Benchmark_comms`, `Benchmark_memory_bandwidth`, … used to qualify a new machine
+- `examples/` — small, readable programmes (`Example_plaquette.cc`, `Example_Mobius_spectrum.cc`) that are the best starting point for learning the API
+
+Each of these directories auto-builds every top-level `.cc` as its own binary via `scripts/filelist`.
+
+Every programme is wrapped in `Grid_init(&argc, &argv)` / `Grid_finalize()` (`Grid/util/Init.h`).

 ## Key Conventions

 - **C++17** is required throughout.
- Template structure: most classes are templated on `<_FImpl>` (fermion impl) or `<Gimpl>` (gauge impl), which encode the representation and precision. Instantiation is controlled by `--enable-fermion-instantiations`.
+- Template structure: most classes are templated on `<_FImpl>` (fermion impl) or `<Gimpl>` (gauge impl), which encode the representation and precision. Instantiation is controlled by `--enable-fermion-instantiations`; the explicit instantiation `.cc` files live under `Grid/qcd/action/fermion/instantiation/<Impl>/` and are selected by `scripts/filelist`.
 - The `RealD`/`RealF`/`ComplexD`/`ComplexF` typedefs are used everywhere; avoid raw `double`/`float`.
- Logging uses `Grid_log`, `Grid_error` macros (from `Grid/log/`); performance-critical paths use the `GRID_TRACE` / timer macros from `Grid/perfmon/`.
+- Use `GRID_ASSERT(cond)` (defined in `Grid/GridStd.h`), not bare `assert` — it prints a Grid-formatted message and aborts cleanly under MPI.
+- Logging is stream-based, not macro-based: `std::cout << GridLogMessage << ... << std::endl;`. Channels declared in `Grid/log/Log.h` include `GridLogError`, `GridLogWarning`, `GridLogDebug`, `GridLogPerformance`, `GridLogIterative`, `GridLogSolver`, `GridLogHMC`, `GridLogComms`, `GridLogMemory`, `GridLogDslash`, `GridLogIRL`, `GridLogMG`. A subset is switched on at runtime with e.g. `--log Error,Warning,Message,Performance,Iterative,Integrator,Debug,Colours` (names given without the `GridLog` prefix).
+- Performance-critical paths use `GRID_TRACE(name)` from `Grid/perfmon/Tracing.h` (compiled out unless `--enable-tracing` selects a backend) and the `GridStopWatch` timers in `Grid/perfmon/Timer.h`.
 - Reductions across MPI ranks go through `GridBase::GlobalSum` / `GlobalMax`; never reduce with bare MPI calls inside library code.
-
-## Skills
-
-`skills/` contains seven user-invocable Claude Code skills encoding deep domain knowledge for HPC work in this repo. Invoke with `/skill-name` or ask Claude to use them by name:
-
-| Skill | When to use |
-|-------|-------------|
-| `gpu-memory-performance` | Bandwidth/occupancy problems — `acceleratorThreads()` pitfalls, `coalescedRead/Write`, fused vs staged HBM patterns |
-| `gpu-runtime-correctness` | Wrong answers, non-deterministic results, premature `q.wait()` returns |
-| `communication-overlap` | Designing GPU+MPI overlap pipelines; replacing broken accelerator-aware MPI paths with host-staged 7-phase pipeline |
-| `mpi-heterogeneous` | Collective hangs, buffer aliasing in `MPI_Sendrecv`, heterogeneous topology bugs |
-| `hang-diagnosis` | Distinguishing ioctl hangs, infinite poll loops, collective deadlocks, and silent wrong-answer races |
-| `correctness-verification` | Reproducibility checksums, double-wait testing, bisecting non-deterministic failures |
-| `compiler-validation` | Confirming compiler/optimisation flags are safe before production runs |
+- Everything lives in `NAMESPACE_BEGIN(Grid)` / `NAMESPACE_END(Grid)` macros; follow the surrounding file rather than writing `namespace Grid { }`.
+- British spelling is used in identifiers and comments (`colour`, `neighbour`, `serialisation`).
@@ -288,21 +288,17 @@ public:
      
      //----------------------------------------------------------------------
      if ( Nm>Nk ) {
-        Glog <<" #Ritz values of poly(A) before shift ["<<Nm<<" values]:"<<std::endl;
-        for(int i=0; i<Nm; ++i){
-          std::cout.precision(8);
-          std::cout << "[" << std::setw(4)<< std::setiosflags(std::ios_base::right) <<i<<"] ";
-          std::cout << "Rval = "<<std::setw(16)<< std::setiosflags(std::ios_base::left)<< eval2[i] << std::endl;
-        }
-
+	//        Glog <<" #Apply shifted QR transformations "<<std::endl;
+        //int k2 = Nk+Nu;
        int k2 = Nk;
-
+      
        Eigen::MatrixXcd BTDM = Eigen::MatrixXcd::Identity(Nm,Nm);
        Q = Eigen::MatrixXcd::Identity(Nm,Nm);
-
+        
        unpackHermitBlockTriDiagMatToEigen(lmd,lme,Nu,Nblock_m,Nm,Nm,BTDM);

        for(int ip=Nk; ip<Nm; ++ip){
+	  Glog << " ip "<<ip<<" / "<<Nm<<std::endl;
          shiftedQRDecompEigen(BTDM,Nu,Nm,eval2[ip],Q);
        }
        
@@ -330,11 +326,11 @@ public:
        Qt = Eigen::MatrixXcd::Identity(Nm,Nm);
        diagonalize(eval2,lmd2,lme2,Nu,Nk,Nm,Qt,grid);
        _sort.push(eval2,Nk);
-        Glog << "#Ritz values of poly(A) after shift ["<<Nk<<" values]:"<<std::endl;
+	//	Glog << "#Ritz value after shift: "<< std::endl;
        for(int i=0; i<Nk; ++i){
-          std::cout.precision(8);
-          std::cout << "[" << std::setw(4)<< std::setiosflags(std::ios_base::right) <<i<<"] ";
-          std::cout << "Rval = "<<std::setw(16)<< std::setiosflags(std::ios_base::left)<< eval2[i] << std::endl;
+	  //	  std::cout.precision(13);
+	  //	  std::cout << "[" << std::setw(4)<< std::setiosflags(std::ios_base::right) <<i<<"] ";
+	  //	  std::cout << "Rval = "<<std::setw(20)<< std::setiosflags(std::ios_base::left)<< eval2[i] << std::endl;
        }
      }
      //----------------------------------------------------------------------
@@ -714,17 +710,11 @@ private:
    //lme[0][L] = beta;
    
    for (int u=0; u<Nu; ++u) {
+      //      Glog << "norm2(w[" << u << "])= "<< norm2(w[u]) << std::endl;
      GRID_ASSERT (!isnan(norm2(w[u])));
-    }
-    // Diagnostic: print alpha (lmd) and beta (lme) block diagonals for this step
-    {
-      Glog << " blk b="<<std::setw(3)<<b<<"  alpha:";
-      for (int u=0; u<Nu; ++u)
-        std::cout << " " << std::setw(10) << std::setprecision(6) << real(lmd[u][L+u]);
-      std::cout << "   |beta|:";
-      for (int u=0; u<Nu; ++u)
-        std::cout << " " << std::setw(10) << std::setprecision(6) << abs(lme[u][L+u]);
-      std::cout << std::endl;
+      for (int k=L+u; k<R; ++k) {
+	//        Glog <<" In block "<< b << "," <<" beta[" << u << "," << k-L << "] = " << lme[u][k] << std::endl;
+      }
    }
    //    Glog << "LinAlg done "<< std::endl;

@@ -63,12 +63,10 @@ void MemoryManager::PrintBytes(void)
  std::cout << " MemoryManager : "<<(total_device>>20)<<" accelerator Mbytes "<<std::endl;
  std::cout << " MemoryManager : "<<(total_host>>20)  <<" cpu         Mbytes "<<std::endl;
  uint64_t cacheBytes;
-  cacheBytes = CacheBytes[Cpu];
-  std::cout << " MemoryManager : "<<(cacheBytes>>20) <<" cpu cache Mbytes "<<std::endl;
-  cacheBytes = CacheBytes[Acc];
-  std::cout << " MemoryManager : "<<(cacheBytes>>20) <<" acc cache Mbytes "<<std::endl;
-  cacheBytes = CacheBytes[Shared];
-  std::cout << " MemoryManager : "<<(cacheBytes>>20) <<" shared cache Mbytes "<<std::endl;
+  cacheBytes = HostCacheBytes();
+  std::cout << " MemoryManager : "<<(cacheBytes>>20) <<" cpu alloc cache Mbytes "<<std::endl;
+  cacheBytes = DeviceCacheBytes();
+  std::cout << " MemoryManager : "<<(cacheBytes>>20) <<" acc alloc cache Mbytes "<<std::endl;
  
 #ifdef GRID_CUDA
  cuda_mem();
@@ -238,6 +238,16 @@ public:
  }
  void AllToAll(int dim  ,void *in,void *out,uint64_t words,uint64_t bytes);
  void AllToAll(void  *in,void *out,uint64_t words         ,uint64_t bytes);
+
+  ////////////////////////////////////////////////////////////////////////////
+  // Variable count all to all. Counts and displacements are in units of
+  // "bytes" sized words and are indexed by rank within this communicator.
+  // For exchanges that are a permutation but do not divide evenly between
+  // ranks; AllToAll above is the uniform count special case.
+  ////////////////////////////////////////////////////////////////////////////
+  void AllToAllV(void *in ,const std::vector<int> &sendcounts,const std::vector<int> &senddispls,
+                 void *out,const std::vector<int> &recvcounts,const std::vector<int> &recvdispls,
+                 uint64_t bytes);
  
  template<class obj> void Broadcast(int root,obj &data)
  {
@@ -945,5 +945,25 @@ void CartesianCommunicator::AllToAll(void  *in,void *out,uint64_t words,uint64_t
  MPI_Alltoall(in,iwords,object,out,iwords,object,communicator);
  MPI_Type_free(&object);
 }
+void CartesianCommunicator::AllToAllV(void *in ,const std::vector<int> &sendcounts,const std::vector<int> &senddispls,
+                                      void *out,const std::vector<int> &recvcounts,const std::vector<int> &recvdispls,
+                                      uint64_t bytes)
+{
+  FlightRecorder::StepLog("AllToAllV");
+  GRID_ASSERT(sendcounts.size()==(size_t)_Nprocessors);
+  GRID_ASSERT(senddispls.size()==(size_t)_Nprocessors);
+  GRID_ASSERT(recvcounts.size()==(size_t)_Nprocessors);
+  GRID_ASSERT(recvdispls.size()==(size_t)_Nprocessors);
+  // MPI counts are "int"; the caller sizes the word to keep them in range
+  int ibytes = bytes;
+  GRID_ASSERT(bytes == (uint64_t)ibytes);
+  MPI_Datatype object;
+  MPI_Type_contiguous(ibytes,MPI_BYTE,&object);
+  MPI_Type_commit(&object);
+  int ierr = MPI_Alltoallv(in ,(int *)&sendcounts[0],(int *)&senddispls[0],object,
+                           out,(int *)&recvcounts[0],(int *)&recvdispls[0],object,communicator);
+  GRID_ASSERT(ierr==0);
+  MPI_Type_free(&object);
+}

 NAMESPACE_END(Grid);
@@ -113,6 +113,17 @@ void CartesianCommunicator::AllToAll(void  *in,void *out,uint64_t words,uint64_t
 {
  bcopy(in,out,bytes*words);
 }
+void CartesianCommunicator::AllToAllV(void *in ,const std::vector<int> &sendcounts,const std::vector<int> &senddispls,
+                                      void *out,const std::vector<int> &recvcounts,const std::vector<int> &recvdispls,
+                                      uint64_t bytes)
+{
+  // Single rank: the exchange degenerates to a copy of our own segment
+  GRID_ASSERT(sendcounts.size()==1);
+  GRID_ASSERT(recvcounts.size()==1);
+  GRID_ASSERT(sendcounts[0]==recvcounts[0]);
+  bcopy((char *)in +(uint64_t)senddispls[0]*bytes,
+        (char *)out+(uint64_t)recvdispls[0]*bytes,bytes*(uint64_t)sendcounts[0]);
+}

 int  CartesianCommunicator::RankWorld(void){return 0;}
 void CartesianCommunicator::Barrier(void){}
@@ -289,7 +289,7 @@ public:
  ///////////////////////////////////////////
  // move constructor
  ///////////////////////////////////////////
-  Lattice(Lattice && r){ 
+  Lattice(Lattice && r) noexcept { 
    this->_grid = r.Grid();
    this->_odata      = r._odata;
    this->_odata_size = r._odata_size;
@@ -330,7 +330,7 @@ public:
  ///////////////////////////////////////////
  // Move assignment possible if same type
  ///////////////////////////////////////////
-  inline Lattice<vobj> & operator = (Lattice<vobj> && r){
+  inline Lattice<vobj> & operator = (Lattice<vobj> && r) noexcept {

    resize(0); // deletes if appropriate
    this->_grid       = r.Grid();
@@ -438,5 +438,11 @@ inline typename vobj::scalar_object sum_gpu_large(const vobj *lat, Integer osite
  result = sumD_gpu_large(lat,osites);
  return result;
 }
+template<class Word> Word checksum_gpu(Word *vec,uint64_t L)
+{
+  Word w;
+  bzero(&w,sizeof(w));
+  return w;
+}

 NAMESPACE_END(Grid);
@@ -2,3 +2,7 @@

 int                    Grid::BinaryIO::latticeWriteMaxRetry = -1;
 Grid::BinaryIO::IoPerf Grid::BinaryIO::lastPerf;
+
+// Target size of a single contiguous file extent under BINARYIO_AGGREGATE.
+// 4MB is around the knee for Lustre; exposed so it can be swept at runtime.
+uint64_t               Grid::BinaryIO::aggregateTargetBytes = 4*1024*1024;
@@ -39,6 +39,7 @@
 #endif

 #include <arpa/inet.h>
+#include <sys/stat.h>
 #include <algorithm>

 NAMESPACE_BEGIN(Grid);
@@ -87,6 +88,7 @@ class BinaryIO {

  static IoPerf lastPerf;
  static int latticeWriteMaxRetry;
+  static uint64_t aggregateTargetBytes;

  /////////////////////////////////////////////////////////////////////////////
  // more byte manipulation helpers
@@ -253,12 +255,392 @@ class BinaryIO {
  // Read or Write distributed lexico array of ANY object to a specific location in file 
  //////////////////////////////////////////////////////////////////////////////////////

+  static const int BINARYIO_AGGREGATE     = 0x20;
  static const int BINARYIO_MASTER_APPEND = 0x10;
  static const int BINARYIO_UNORDERED     = 0x08;
  static const int BINARYIO_LEXICOGRAPHIC = 0x04;
  static const int BINARYIO_READ          = 0x02;
  static const int BINARYIO_WRITE         = 0x01;

+#ifdef USE_MPI_IO
+  /////////////////////////////////////////////////////////////////////////////
+  // Aggregation: self controlled transposition onto an I/O friendly layout.
+  //
+  // Under BINARYIO_LEXICOGRAPHIC the subarray file view handed to MPI-IO has
+  // contiguous runs of only lLattice[0]*sizeof(fobj) bytes -- a few KB for
+  // typical local volumes. Rather than rely on collective buffering to repair
+  // that, redistribute the payload ourselves so every rank owns a contiguous
+  // range of the global lexicographic site ordering, then issue large plain
+  // contiguous writes.
+  //
+  // "Un-splitting" the nunsplit fastest dimensions means the row of ranks
+  // sharing the remaining process coordinates collectively owns whole global
+  // hyperplanes. All data movement is then confined to that row communicator.
+  // Every rank still owns exactly lSites() sites afterwards, so the exchange is
+  // a pure permutation and needs no divisibility condition on the process grid.
+  /////////////////////////////////////////////////////////////////////////////
+  struct AggregationPlan {
+    int      nunsplit{0};  // number of fastest dimensions un-split
+    int      rowsize{0};   // ranks in the aggregation (row) communicator
+    int      rowrank{0};   // our logical (lexicographic) index within the row
+    uint64_t lsites{0};    // sites per rank -- invariant under the permutation
+    uint64_t chunk{0};     // sites in one globally contiguous run owned by the row
+    std::unique_ptr<CartesianCommunicator> rowcomm;
+    // counts and displacements are indexed by rank within rowcomm
+    std::vector<int>      sendcounts, senddispls, recvcounts, recvdispls;
+    std::vector<uint64_t> scatter;     // recv slot -> slot in the aggregated buffer
+    std::vector<uint64_t> extentGsite; // global lex site index of extent start
+    std::vector<uint64_t> extentLocal; // offset of extent within aggregated buffer
+    std::vector<uint64_t> extentSites; // sites in this extent
+  };
+
+  static inline void BuildAggregationPlan(GridBase *grid,uint64_t fobjSize,AggregationPlan &p)
+  {
+    int ndim           = grid->Dimensions();
+    Coordinate psizes  = grid->ProcessorGrid();
+    Coordinate pcoor   = grid->ThisProcessorCoor();
+    Coordinate gLattice= grid->GlobalDimensions();
+    Coordinate lLattice= grid->LocalDimensions();
+    Coordinate lstart  = grid->LocalStarts();
+
+    uint64_t lsites = grid->lSites();
+    p.lsites = lsites;
+
+    //////////////////////////////////////////////////////////////////////////
+    // Un-splitting dims 0..k-1 gives the row a contiguous run of
+    //     chunk(k) = prod_{d<k} gLattice[d] * lLattice[k]
+    // sites, and each rank writes extents of min(chunk,lsites). Take the
+    // smallest k that reaches the target so we disturb as few dimensions --
+    // and move as little data -- as possible.
+    //////////////////////////////////////////////////////////////////////////
+    int k = ndim-1;
+    for(int trial=1; trial<ndim; trial++){
+      uint64_t chunk = lLattice[trial];
+      for(int d=0; d<trial; d++) chunk *= gLattice[d];
+      if ( std::min(chunk,lsites)*fobjSize >= aggregateTargetBytes ) { k = trial; break; }
+    }
+    p.nunsplit = k;
+
+    //////////////////////////////////////////////////////////////////////////
+    // The box the row collectively owns, expressed in global coordinates.
+    // Restricting the global lexicographic order to this box preserves the
+    // ordering, so the row index below is monotone in the global index.
+    //////////////////////////////////////////////////////////////////////////
+    Coordinate B(ndim), S(ndim);
+    for(int d=0; d<ndim; d++){
+      if ( d<k ) { B[d] = gLattice[d]; S[d] = 0;         }
+      else       { B[d] = lLattice[d]; S[d] = lstart[d]; }
+    }
+
+    uint64_t chunk = lLattice[k];
+    for(int d=0; d<k; d++) chunk *= gLattice[d];
+    p.chunk = chunk;
+
+    //////////////////////////////////////////////////////////////////////////
+    // Row communicator: the ranks sharing the process coordinates of the slow
+    // (still split) dimensions. This is the sub-division the Cartesian
+    // communicator already performs for AllToAll(dim,...), widened from one
+    // dimension to the k fastest.
+    //////////////////////////////////////////////////////////////////////////
+    Coordinate row(ndim,1);
+    for(int d=0; d<k; d++) row[d] = psizes[d];
+    int srank;
+    p.rowcomm.reset(new CartesianCommunicator(row,*grid,srank));
+    p.rowsize = p.rowcomm->ProcessorCount();
+
+    //////////////////////////////////////////////////////////////////////////
+    // Our logical index in the row is the forward lexicographic index of the
+    // un-split process coordinates, so that increasing logical index means
+    // increasing global lexicographic position in the file. The communicator
+    // numbers its own ranks by the reversed (MPI) convention, so build the map
+    // between the two rather than assuming either.
+    //////////////////////////////////////////////////////////////////////////
+    int64_t logical=0, lstride=1;
+    for(int d=0; d<k; d++){ logical += pcoor[d]*lstride; lstride *= psizes[d]; }
+    GRID_ASSERT(lstride == (int64_t)p.rowsize);
+    p.rowrank = (int)logical;
+
+    std::vector<uint64_t> commOf(p.rowsize,0);
+    commOf[p.rowrank] = (uint64_t)p.rowcomm->ThisRank();
+    p.rowcomm->GlobalSumVector(&commOf[0],p.rowsize);
+
+    uint64_t mystart = (uint64_t)p.rowrank * lsites;
+    uint64_t myend   = mystart + lsites;
+
+    Coordinate lcoor(ndim), bcoor(ndim), gcoor(ndim);
+
+    //////////////////////////////////////////////////////////////////////////
+    // Send side. Walking our local sites in local lexicographic order walks
+    // the row index monotonically, so the send buffer is iodata untouched and
+    // we need only the per destination counts.
+    //////////////////////////////////////////////////////////////////////////
+    std::vector<int> sendLogical(p.rowsize,0);
+    for(uint64_t L=0; L<lsites; L++){
+      Lexicographic::CoorFromIndex(lcoor,L,lLattice);
+      for(int d=0; d<ndim; d++) bcoor[d] = (d<k) ? (lstart[d]+lcoor[d]) : lcoor[d];
+      int64_t ri; Lexicographic::IndexFromCoor(bcoor,ri,B);
+      sendLogical[ ri/(int64_t)lsites ]++;
+    }
+    p.sendcounts.assign(p.rowsize,0);
+    p.senddispls.assign(p.rowsize,0);
+    { int64_t disp=0;
+      for(int d=0; d<p.rowsize; d++){          // send buffer is in logical order
+        int c = (int)commOf[d];
+        p.sendcounts[c] = sendLogical[d];
+        p.senddispls[c] = (int)disp;
+        disp += sendLogical[d];
+      }
+    }
+
+    //////////////////////////////////////////////////////////////////////////
+    // Receive side. For each slot of our aggregated range work out which rank
+    // of the row owns it. Within one source the slots arrive in increasing row
+    // index order, which is the order the source sends them in.
+    //////////////////////////////////////////////////////////////////////////
+    std::vector<int> recvLogical(p.rowsize,0), recvDisplLogical(p.rowsize,0);
+    std::vector<int> source(lsites);
+    for(uint64_t pos=0; pos<lsites; pos++){
+      Lexicographic::CoorFromIndex(bcoor,(int64_t)(mystart+pos),B);
+      int64_t j=0, jstride=1;
+      for(int d=0; d<k; d++){ j += (bcoor[d]/lLattice[d])*jstride; jstride *= psizes[d]; }
+      source[pos] = (int)j;
+      recvLogical[j]++;
+    }
+    p.recvcounts.assign(p.rowsize,0);
+    p.recvdispls.assign(p.rowsize,0);
+    { int64_t disp=0;
+      for(int s=0; s<p.rowsize; s++){         // recv buffer is in logical order
+        int c = (int)commOf[s];
+        recvDisplLogical[s] = (int)disp;
+        p.recvcounts[c] = recvLogical[s];
+        p.recvdispls[c] = (int)disp;
+        disp += recvLogical[s];
+      }
+    }
+    p.scatter.resize(lsites);
+    {
+      std::vector<int> fill(p.rowsize,0);
+      for(uint64_t pos=0; pos<lsites; pos++){
+        int j = source[pos];
+        p.scatter[ recvDisplLogical[j] + fill[j]++ ] = pos;
+      }
+    }
+
+    //////////////////////////////////////////////////////////////////////////
+    // The two sides are derived independently; make them check each other.
+    //////////////////////////////////////////////////////////////////////////
+    {
+      std::vector<uint64_t> sendc(p.rowsize),recvc(p.rowsize);
+      for(int c=0;c<p.rowsize;c++) sendc[c]=(uint64_t)p.sendcounts[c];
+      p.rowcomm->AllToAll(&sendc[0],&recvc[0],1,sizeof(uint64_t));
+      for(int c=0;c<p.rowsize;c++) GRID_ASSERT((int)recvc[c]==p.recvcounts[c]);
+    }
+
+    //////////////////////////////////////////////////////////////////////////
+    // Decompose our range into globally contiguous file extents.
+    //////////////////////////////////////////////////////////////////////////
+    for(uint64_t c = mystart/chunk; c <= (myend-1)/chunk; c++){
+      uint64_t lo = std::max(mystart, c*chunk);
+      uint64_t hi = std::min(myend, (c+1)*chunk);
+      Lexicographic::CoorFromIndex(bcoor,(int64_t)(c*chunk),B);
+      for(int d=0;d<ndim;d++) gcoor[d] = (d<k) ? bcoor[d] : bcoor[d]+S[d];
+      int64_t gbase; Lexicographic::IndexFromCoor(gcoor,gbase,gLattice);
+      p.extentGsite.push_back( (uint64_t)gbase + (lo - c*chunk) );
+      p.extentLocal.push_back( lo - mystart );
+      p.extentSites.push_back( hi - lo );
+    }
+  }
+
+  static inline void ReportAggregationPlan(GridBase *grid,const AggregationPlan &p,uint64_t fobjSize,const char *what)
+  {
+    if ( !grid->IsBoss() ) return;
+    std::cout << GridLogMessage << "IOobject: aggregate " << what
+              << " un-splitting " << p.nunsplit << " fastest dimensions, row of "
+              << p.rowsize << " ranks" << std::endl;
+    std::cout << GridLogMessage << "IOobject: aggregate " << p.extentSites.size()
+              << " extent(s)/rank, first " << p.extentSites[0]*fobjSize/1024./1024. << " MB"
+              << " (target " << aggregateTargetBytes/1024./1024. << " MB)" << std::endl;
+    std::cout << GridLogMessage << "IOobject: aggregate buffer overhead "
+              << p.lsites*fobjSize/1024./1024. << " MB/rank" << std::endl;
+  }
+
+  ////////////////////////////////////////////////////////////////////////////
+  // Stage timings.  The interesting quantity is the slowest rank, since every
+  // stage is followed sooner or later by a synchronisation, so reduce with
+  // GlobalMax rather than reporting whatever the boss happened to see.
+  ////////////////////////////////////////////////////////////////////////////
+  static inline void ReportStages(GridBase *grid,const char *what,
+                                  const std::vector<const char *> &names,
+                                  std::vector<RealD> &useconds)
+  {
+    GRID_ASSERT(names.size()==useconds.size());
+    for(uint64_t i=0;i<useconds.size();i++) grid->GlobalMax(useconds[i]);
+    if ( grid->IsBoss() ) {
+      std::cout << GridLogMessage << "IOobject: aggregate " << what << " stages (max over ranks, s):";
+      for(uint64_t i=0;i<names.size();i++)
+        std::cout << " " << names[i] << " " << useconds[i]/1.0e6;
+      std::cout << std::endl;
+    }
+  }
+
+  template<class fobj>
+  static inline void AggregateExchange(GridBase *grid,AggregationPlan &p,std::vector<fobj> &iodata,
+                                       std::vector<fobj> &aggregated,int forward)
+  {
+    uint64_t lsites = p.lsites;
+    GridStopWatch talloc,tperm,tcomm;
+
+    talloc.Start();
+    std::vector<fobj> tmp(lsites);
+    talloc.Stop();
+
+    if ( forward ) { // iodata (local order) -> aggregated (lexicographic order)
+      tcomm.Start();
+      p.rowcomm->AllToAllV(&iodata[0],p.sendcounts,p.senddispls,
+                           &tmp[0],   p.recvcounts,p.recvdispls,sizeof(fobj));
+      tcomm.Stop();
+      tperm.Start();
+      thread_for(s,lsites,{ aggregated[p.scatter[s]] = tmp[s]; });
+      tperm.Stop();
+    } else {         // aggregated -> iodata, the exact mirror
+      tperm.Start();
+      thread_for(s,lsites,{ tmp[s] = aggregated[p.scatter[s]]; });
+      tperm.Stop();
+      tcomm.Start();
+      p.rowcomm->AllToAllV(&tmp[0],   p.recvcounts,p.recvdispls,
+                           &iodata[0],p.sendcounts,p.senddispls,sizeof(fobj));
+      tcomm.Stop();
+    }
+
+    std::vector<RealD> us = { (RealD)talloc.useconds(), (RealD)tperm.useconds(), (RealD)tcomm.useconds() };
+    ReportStages(grid,forward?"exchange (write)":"exchange (read)",
+                 {"alloc","permute","alltoallv"},us);
+  }
+
+  template<class fobj>
+  static inline void AggregateWrite(GridBase *grid,AggregationPlan &p,std::vector<fobj> &aggregated,
+                                    std::string file,uint64_t offset)
+  {
+    //////////////////////////////////////////////////////////////////////////
+    // All ranks write concurrently into a shared file, so the file must exist
+    // before any of them open it for update, but it does NOT have to be the
+    // right length first: the extents tile the record exactly, so writing them
+    // extends a short file to precisely offset+payload.
+    //
+    // Records are created in sequence, so this payload ends the file: the
+    // length must end up precisely offset+payload. Anything beyond is left
+    // over from whatever the file previously held and must not survive -- a
+    // shorter new record written over a longer old one would otherwise leave
+    // a trailing fragment of the previous contents masquerading as data.
+    // That is the only case needing a truncate, so stat first and truncate
+    // afterwards only when the size actually came out wrong. Measured on
+    // Frontier, an unconditional truncate up front cost 0.22 to 5.4 s per
+    // record -- 15 to 25% of a 19 GB write and 100% of a small one -- while
+    // create, open and close together cost a few milliseconds. It is per
+    // record, so multi record files do not amortise it away.
+    //
+    // ::truncate is used because the C++ standard library cannot express this.
+    // std::filebuf has no length operation at all; ios::trunc only truncates to
+    // zero; seeking past the end and writing a byte can grow a file but never
+    // shrink one; and there is no portable way to recover a descriptor from a
+    // stream in order to call ftruncate. C++17 does finally offer
+    // std::filesystem::resize_file, but that would be Grid's first <filesystem>
+    // dependency and needs -lstdc++fs on the older toolchains still in use.
+    //////////////////////////////////////////////////////////////////////////
+    GridStopWatch tcreate,ttrunc,tbar,topen,twrite,tclose,tskew;
+    uint64_t need = offset + (uint64_t)grid->_gsites*sizeof(fobj);
+
+    tcreate.Start();
+    if ( grid->IsBoss() ) {
+      // opening for update needs the file to exist; create one only if not
+      std::fstream probe(file,std::ios::binary|std::ios::out|std::ios::in);
+      if ( !probe.is_open() ) {
+        std::ofstream create(file,std::ios::binary|std::ios::out);
+        create.close();
+      }
+    }
+    tcreate.Stop();
+
+    tbar.Start();
+    grid->Barrier();
+    tbar.Stop();
+
+    std::ofstream fout;
+    fout.exceptions( std::fstream::failbit | std::fstream::badbit );
+    try {
+      topen.Start();
+      fout.open(file,std::ios::binary|std::ios::out|std::ios::in);
+      topen.Stop();
+      twrite.Start();
+      for(uint64_t e=0;e<p.extentSites.size();e++){
+        fout.seekp(offset + p.extentGsite[e]*sizeof(fobj));
+        fout.write((char *)&aggregated[p.extentLocal[e]],p.extentSites[e]*sizeof(fobj));
+      }
+      twrite.Stop();
+      tclose.Start();
+      fout.close();   // flushes the stream buffer; does not force writeback
+      tclose.Stop();
+    } catch (const std::fstream::failure& exc) {
+      std::cout << GridLogError << "Error in aggregate write to " << file << std::endl;
+      std::cout << GridLogError << "Exception description: " << exc.what() << std::endl;
+      GridAbort();
+    }
+
+    ////////////////////////////////////////////////////////////////////////
+    // Timed apart from the truncate that follows it.  seek+write above is the
+    // slowest rank; this barrier is what the fastest rank then waits, so the
+    // pair separates the write cost from the spread across ranks. Folding it
+    // into the truncate makes a millisecond stat look like a second.
+    ////////////////////////////////////////////////////////////////////////
+    tskew.Start();
+    grid->Barrier();                 // every extent must be on its way first
+    tskew.Stop();
+
+    ttrunc.Start();
+    if ( grid->IsBoss() ) {
+      struct stat sb;
+      int ierr = ::stat(file.c_str(),&sb);
+      GRID_ASSERT(ierr==0);
+      if ( (uint64_t)sb.st_size != need ) {   // only when a longer record preceded us
+        ierr = ::truncate(file.c_str(),(off_t)need);
+        GRID_ASSERT(ierr==0);
+      }
+    }
+    grid->Barrier();
+    ttrunc.Stop();
+
+    std::vector<RealD> us = { (RealD)tcreate.useconds(), (RealD)tbar.useconds(),
+                              (RealD)topen.useconds(),   (RealD)twrite.useconds(),
+                              (RealD)tclose.useconds(),  (RealD)tskew.useconds(),
+                              (RealD)ttrunc.useconds() };
+    ReportStages(grid,"write",{"create","barrier","open","seek+write","close","skew","stat+truncate"},us);
+  }
+
+  template<class fobj>
+  static inline void AggregateRead(GridBase *grid,AggregationPlan &p,std::vector<fobj> &aggregated,
+                                   std::string file,uint64_t offset)
+  {
+    GridStopWatch topen,tread,tclose;
+    std::ifstream fin;
+    topen.Start();
+    fin.open(file,std::ios::binary|std::ios::in);
+    topen.Stop();
+    tread.Start();
+    for(uint64_t e=0;e<p.extentSites.size();e++){
+      fin.seekg(offset + p.extentGsite[e]*sizeof(fobj));
+      fin.read((char *)&aggregated[p.extentLocal[e]],p.extentSites[e]*sizeof(fobj));
+      GRID_ASSERT(fin.fail()==0);
+    }
+    tread.Stop();
+    tclose.Start();
+    fin.close();
+    tclose.Stop();
+
+    std::vector<RealD> us = { (RealD)topen.useconds(), (RealD)tread.useconds(), (RealD)tclose.useconds() };
+    ReportStages(grid,"read",{"open","seek+read","close"},us);
+  }
+#endif
+
  template<class word,class fobj>
  static inline void IOobject(word w,
 			      GridBase *grid,
@@ -302,6 +684,18 @@ class BinaryIO {
      lStart[d] = 0;
    }

+    //////////////////////////////////////////////////////////////////////////////
+    // Aggregate the lexicographic layout onto contiguous per rank extents
+    // ourselves rather than leaving it to MPI-IO collective buffering
+    //////////////////////////////////////////////////////////////////////////////
+    int aggregate = (control & BINARYIO_AGGREGATE)
+                 && (control & BINARYIO_LEXICOGRAPHIC)
+                 && !(control & BINARYIO_MASTER_APPEND)
+                 && (nrank > 1);
+#ifndef USE_MPI_IO
+    GRID_ASSERT(aggregate==0); // BINARYIO_AGGREGATE requires MPI
+#endif
+
 #ifdef USE_MPI_IO
    std::vector<int> distribs(ndim,MPI_DISTRIBUTE_BLOCK);
    std::vector<int> dargs   (ndim,MPI_DISTRIBUTE_DFLT_DARG);
@@ -329,6 +723,8 @@ class BinaryIO {
    ierr = MPI_Type_contiguous(numword,mpiword,&mpiObject);    GRID_ASSERT(ierr==0);
    ierr = MPI_Type_commit(&mpiObject);

+    // The subarray view is what aggregation exists to avoid; do not build it
+    if ( !aggregate ) {
    //////////////////////////////////////////////////////////////////////////////
    // File global array data type
    //////////////////////////////////////////////////////////////////////////////
@@ -340,6 +736,7 @@ class BinaryIO {
    //////////////////////////////////////////////////////////////////////////////
    ierr=MPI_Type_create_subarray(ndim,&lLattice[0],&lLattice[0],&lStart[0],MPI_ORDER_FORTRAN, mpiObject,&localArray);    GRID_ASSERT(ierr==0);
    ierr=MPI_Type_commit(&localArray);    GRID_ASSERT(ierr==0);
+    }
 #endif

    //////////////////////////////////////////////////////////////////////////////
@@ -358,7 +755,19 @@ class BinaryIO {

      timer.Start();

-      if ( (control & BINARYIO_LEXICOGRAPHIC) && (nrank > 1) ) {
+      if ( aggregate ) {
+#ifdef USE_MPI_IO
+        std::cout<< GridLogMessage<<"IOobject: aggregate read I/O "<< file<< std::endl;
+        AggregationPlan plan;
+        BuildAggregationPlan(grid,sizeof(fobj),plan);
+        ReportAggregationPlan(grid,plan,sizeof(fobj),"read");
+        std::vector<fobj> aggregated(lsites);
+        AggregateRead(grid,plan,aggregated,file,offset);
+        AggregateExchange(grid,plan,iodata,aggregated,0);
+#else
+        GRID_ASSERT(0);
+#endif
+      } else if ( (control & BINARYIO_LEXICOGRAPHIC) && (nrank > 1) ) {
 #ifdef USE_MPI_IO
 	std::cout<< GridLogMessage<<"IOobject: MPI read I/O "<< file<< std::endl;
 	ierr=MPI_File_open(grid->communicator,(char *) file.c_str(), MPI_MODE_RDONLY, MPI_INFO_NULL, &fh);    GRID_ASSERT(ierr==0);
@@ -416,7 +825,25 @@ class BinaryIO {
      grid->Barrier();

      timer.Start();
-      if ( (control & BINARYIO_LEXICOGRAPHIC) && (nrank > 1) ) {
+      if ( aggregate ) {
+#ifdef USE_MPI_IO
+        std::cout << GridLogMessage <<"IOobject: aggregate write I/O " << file << std::endl;
+        AggregationPlan plan;
+        BuildAggregationPlan(grid,sizeof(fobj),plan);
+        ReportAggregationPlan(grid,plan,sizeof(fobj),"write");
+        std::vector<fobj> aggregated(lsites);
+        AggregateExchange(grid,plan,iodata,aggregated,1);
+        AggregateWrite(grid,plan,aggregated,file,offset);
+        ////////////////////////////////////////////////////////////////////////
+        // Not every rank ends at the end of the payload, so the position can
+        // not be recovered from a file handle. Callers (Lime record chaining)
+        // rely on this being the first byte past the record.
+        ////////////////////////////////////////////////////////////////////////
+        offset = offset + (uint64_t)grid->_gsites*sizeof(fobj);
+#else
+        GRID_ASSERT(0);
+#endif
+      } else if ( (control & BINARYIO_LEXICOGRAPHIC) && (nrank > 1) ) {
 #ifdef USE_MPI_IO
        std::cout << GridLogMessage <<"IOobject: MPI write I/O " << file << std::endl;
        ierr = MPI_File_open(grid->communicator, (char *)file.c_str(), MPI_MODE_RDWR | MPI_MODE_CREATE, MPI_INFO_NULL, &fh);
@@ -461,12 +888,26 @@ class BinaryIO {
        
 	std::ofstream fout; 
 	fout.exceptions ( std::fstream::failbit | std::fstream::badbit );
+
+        ////////////////////////////////////////////////////////////////////
+        // Grid's model is that the boss rank performs the metadata
+        // operations and every other rank only seeks and writes into a file
+        // that already exists. Opening with ios::out on all ranks broke that:
+        // it is O_TRUNC, so a rank opening late truncated the file back to
+        // zero after an earlier rank had written its segment, leaving a hole
+        // in its place. The barriers around this block are outside it and do
+        // not order the opens against the writes. Let the boss create and
+        // empty the file, then everyone opens for update only. Same resulting
+        // length, one metadata operation instead of one per rank, no race.
+        ////////////////////////////////////////////////////////////////////
+        if ( !offset && grid->IsBoss() ) { // offset zero: this record starts the file
+          std::ofstream create(file,std::ios::binary|std::ios::out);
+          create.close();
+        }
+        grid->Barrier();
+
 	try {
-	  if (offset) { // Must already exist and contain data
-	    fout.open(file,std::ios::binary|std::ios::out|std::ios::in);
-	  } else {     // Allow create
-	    fout.open(file,std::ios::binary|std::ios::out);
-	  }
+	  fout.open(file,std::ios::binary|std::ios::out|std::ios::in);
 	} catch (const std::fstream::failure& exc) {
 	  std::cout << GridLogError << "Error in opening the file " << file << " for output" <<std::endl;
 	  std::cout << GridLogError << "Exception description: " << exc.what() << std::endl;
@@ -477,7 +918,7 @@ class BinaryIO {
 	  exit(1);
 #endif
 	}
-	
+
 	if ( control & BINARYIO_MASTER_APPEND )  {
 	  try {
 	    fout.seekp(0,fout.end);
@@ -548,7 +989,7 @@ class BinaryIO {
 				       uint32_t &nersc_csum,
 				       uint32_t &scidac_csuma,
 				       uint32_t &scidac_csumb,
-				       int control=BINARYIO_LEXICOGRAPHIC
+				       int control=BINARYIO_LEXICOGRAPHIC|BINARYIO_AGGREGATE
 				       )
  {
    typedef typename vobj::scalar_object sobj;
@@ -587,7 +1028,7 @@ class BinaryIO {
 					  uint32_t &nersc_csum,
 					  uint32_t &scidac_csuma,
 					  uint32_t &scidac_csumb,
-					  int control=BINARYIO_LEXICOGRAPHIC)
+					  int control=BINARYIO_LEXICOGRAPHIC|BINARYIO_AGGREGATE)
  {
    typedef typename vobj::scalar_object sobj;
    typedef typename vobj::Realified::scalar_type word;    word w=0;
@@ -674,7 +1115,7 @@ class BinaryIO {
    std::cout << GridLogMessage << "RNG read I/O on file " << file << std::endl;

    std::vector<RNGstate> iodata(lsites);
-    IOobject(w,grid,iodata,file,offset,format,BINARYIO_READ|BINARYIO_LEXICOGRAPHIC,
+    IOobject(w,grid,iodata,file,offset,format,BINARYIO_READ|BINARYIO_LEXICOGRAPHIC|BINARYIO_AGGREGATE,
 	     nersc_csum,scidac_csuma,scidac_csumb);

    timer.Start();
@@ -753,7 +1194,7 @@ class BinaryIO {
    });
    timer.Stop();

-    IOobject(w,grid,iodata,file,offset,format,BINARYIO_WRITE|BINARYIO_LEXICOGRAPHIC,
+    IOobject(w,grid,iodata,file,offset,format,BINARYIO_WRITE|BINARYIO_LEXICOGRAPHIC|BINARYIO_AGGREGATE,
 	     nersc_csum,scidac_csuma,scidac_csumb);
    iodata.resize(1);
    {
@@ -51,8 +51,8 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #endif
 #ifdef __x86_64__
 #ifdef GRID_CUDA
-//accelerator_inline uint64_t __rdtsc(void) {  return 0; }
-//accelerator_inline uint64_t __rdpmc(int ) {  return 0; }
+accelerator_inline uint64_t __rdtsc(void) {  return 0; }
+accelerator_inline uint64_t __rdpmc(int ) {  return 0; }
 #else
 #include <x86intrin.h>
 #endif
@@ -93,8 +93,7 @@ inline uint64_t cyclecount(void){
 }
 #elif defined __x86_64__
 inline uint64_t cyclecount(void){ 
-  uint64_t ret = __rdtsc();
-  return (uint64_t)ret;
+  return (uint64_t)0;
 }
 #else

@@ -115,10 +115,10 @@ accelerator_inline ComplexF adj(const ComplexF& r ){ return(conjugate(r)); }

 #if defined(GRID_CUDA) || defined(GRID_HIP)
 //Provide for convenience
-accelerator_inline std::complex<double> conjugate(const std::complex<double>& r){ return(conj(r)); }
-accelerator_inline std::complex<float>  conjugate(const std::complex<float>& r) { return(conj(r)); }
-accelerator_inline std::complex<double> adj(const std::complex<double>& r)      { return(conj(r)); }
-accelerator_inline std::complex<float>  adj(const std::complex<float>& r)       { return(conj(r)); }
+inline std::complex<double> conjugate(const std::complex<double>& r){ return(conj(r)); }
+inline std::complex<float>  conjugate(const std::complex<float>& r) { return(conj(r)); }
+inline std::complex<double> adj(const std::complex<double>& r)      { return(conj(r)); }
+inline std::complex<float>  adj(const std::complex<float>& r)       { return(conj(r)); }
 #endif

 accelerator_inline RealF real(const RealF  & r){ return r; }
@@ -33,6 +33,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>

 using namespace std;
 using namespace Grid;
+ ;

 int main (int argc, char ** argv)
 {
@@ -96,49 +97,19 @@ int main (int argc, char ** argv)
  RealD c2=-1.0/24.0;
  RealD u0=1.0;
  ImprovedStaggeredFermionD Ds(Umu,Umu,Grid,RBGrid,mass,c1,c2,u0,params);
-  NaiveStaggeredFermionD Dn(Umu,Grid,RBGrid,mass,c1,u0,params);
-
-  FermionField src_e(&RBGrid); pickCheckerboard(Even, src_e, src);
-  FermionField src_o(&RBGrid); pickCheckerboard(Odd,  src_o, src);
-  FermionField res_e(&RBGrid); res_e = Zero();
-
+  
+  std::cout<<GridLogMessage << "Calling Ds"<<std::endl;
  int ncall=1000;
-
-  // ImprovedStaggered Dhop
-  for(int i=0;i<ncall;i++) Ds.Dhop(src,result,0);
  double t0=usecond();
-  for(int i=0;i<ncall;i++) Ds.Dhop(src,result,0);
+  for(int i=0;i<ncall;i++){
+    Ds.Dhop(src,result,0);
+  }
  double t1=usecond();
-  double flops=(16*(3*(6+8+8)) + 15*3*2)*volume*ncall; // 1146 flops/site
-  std::cout<<GridLogMessage << ncall << " ImprovedStaggered Dhop calls in "<< (t1-t0)<<" us"<<std::endl;
-  std::cout<<GridLogMessage << "ImprovedStaggered Dhop   mflop/s =   "<< flops/(t1-t0)<<std::endl;
-
-  // ImprovedStaggered DhopEO
-  for(int i=0;i<ncall;i++) Ds.DhopEO(src_o,res_e,0);
-  t0=usecond();
-  for(int i=0;i<ncall;i++) Ds.DhopEO(src_o,res_e,0);
-  t1=usecond();
-  flops=(16*(3*(6+8+8)) + 15*3*2)*(volume/2)*ncall;
-  std::cout<<GridLogMessage << ncall << " ImprovedStaggered DhopEO calls in "<< (t1-t0)<<" us"<<std::endl;
-  std::cout<<GridLogMessage << "ImprovedStaggered DhopEO mflop/s =   "<< flops/(t1-t0)<<std::endl;
-
-  // NaiveStaggered Dhop
-  for(int i=0;i<ncall;i++) Dn.Dhop(src,result,0);
-  t0=usecond();
-  for(int i=0;i<ncall;i++) Dn.Dhop(src,result,0);
-  t1=usecond();
-  flops=(8*(3*(6+8+8)) + 7*3*2)*volume*ncall;
-  std::cout<<GridLogMessage << ncall << " NaiveStaggered Dhop calls in "<< (t1-t0)<<" us"<<std::endl;
-  std::cout<<GridLogMessage << "NaiveStaggered    Dhop   mflop/s =   "<< flops/(t1-t0)<<std::endl;
-
-  // NaiveStaggered DhopEO
-  for(int i=0;i<ncall;i++) Dn.DhopEO(src_o,res_e,0);
-  t0=usecond();
-  for(int i=0;i<ncall;i++) Dn.DhopEO(src_o,res_e,0);
-  t1=usecond();
-  flops=(8*(3*(6+8+8)) + 7*3*2)*(volume/2)*ncall;
-  std::cout<<GridLogMessage << ncall << " NaiveStaggered DhopEO calls in "<< (t1-t0)<<" us"<<std::endl;
-  std::cout<<GridLogMessage << "NaiveStaggered    DhopEO mflop/s =   "<< flops/(t1-t0)<<std::endl;
+  double flops=(16*(3*(6+8+8)) + 15*3*2)*volume*ncall; // == 66*16 +  == 1146
+  
+  std::cout<<GridLogMessage << "Called Ds"<<std::endl;
+  std::cout<<GridLogMessage << "norm result "<< norm2(result)<<std::endl;
+  std::cout<<GridLogMessage << "mflop/s =   "<< flops/(t1-t0)<<std::endl;

  Grid_finalize();
 }
@@ -96,45 +96,22 @@ int main (int argc, char ** argv)
  RealD c2=-1.0/24.0;
  RealD u0=1.0;
  ImprovedStaggeredFermionF Ds(Umu,Umu,Grid,RBGrid,mass,c1,c2,u0,params);
-  NaiveStaggeredFermionF    Dn(Umu,Grid,RBGrid,mass,c1,u0,params);
-
-  FermionField src_e(&RBGrid); pickCheckerboard(Even, src_e, src);
-  FermionField src_o(&RBGrid); pickCheckerboard(Odd,  src_o, src);
-  FermionField res_e(&RBGrid); res_e = Zero();
-
-  int ncall=10000;
-
-  // ImprovedStaggered Dhop
-  for(int i=0;i<ncall;i++) Ds.Dhop(src,result,0);
+  
+  std::cout<<GridLogMessage << "Calling Ds"<<std::endl;
+  int ncall=1000;
+  for(int i=0;i<ncall;i++){
+    Ds.Dhop(src,result,0);
+  }
  double t0=usecond();
-  for(int i=0;i<ncall;i++) Ds.Dhop(src,result,0);
+  for(int i=0;i<ncall;i++){
+    Ds.Dhop(src,result,0);
+  }
  double t1=usecond();
-  double flops=(16*(3*(6+8+8)) + 15*3*2)*volume*ncall; // 1146 flops/site
-  std::cout<<GridLogMessage << "ImprovedStaggered Dhop   mflop/s =   "<< flops/(t1-t0)<<std::endl;
-
-  // ImprovedStaggered DhopEO
-  for(int i=0;i<ncall;i++) Ds.DhopEO(src_o,res_e,0);
-  t0=usecond();
-  for(int i=0;i<ncall;i++) Ds.DhopEO(src_o,res_e,0);
-  t1=usecond();
-  flops=(16*(3*(6+8+8)) + 15*3*2)*(volume/2)*ncall;
-  std::cout<<GridLogMessage << "ImprovedStaggered DhopEO mflop/s =   "<< flops/(t1-t0)<<std::endl;
-
-  // NaiveStaggered Dhop
-  for(int i=0;i<ncall;i++) Dn.Dhop(src,result,0);
-  t0=usecond();
-  for(int i=0;i<ncall;i++) Dn.Dhop(src,result,0);
-  t1=usecond();
-  flops=(8*(3*(6+8+8)) + 7*3*2)*volume*ncall;
-  std::cout<<GridLogMessage << "NaiveStaggered    Dhop   mflop/s =   "<< flops/(t1-t0)<<std::endl;
-
-  // NaiveStaggered DhopEO
-  for(int i=0;i<ncall;i++) Dn.DhopEO(src_o,res_e,0);
-  t0=usecond();
-  for(int i=0;i<ncall;i++) Dn.DhopEO(src_o,res_e,0);
-  t1=usecond();
-  flops=(8*(3*(6+8+8)) + 7*3*2)*(volume/2)*ncall;
-  std::cout<<GridLogMessage << "NaiveStaggered    DhopEO mflop/s =   "<< flops/(t1-t0)<<std::endl;
+  double flops=(16*(3*(6+8+8)) + 15*3*2)*volume*ncall; // == 66*16 +  == 1146
+  
+  std::cout<<GridLogMessage << "Called Ds"<<std::endl;
+  std::cout<<GridLogMessage << "norm result "<< norm2(result)<<std::endl;
+  std::cout<<GridLogMessage << "mflop/s =   "<< flops/(t1-t0)<<std::endl;

  Grid_finalize();
 }
@@ -716,161 +716,6 @@ public:
    return mflops_best;
  }

-  static double NaiveStaggered(int L)
-  {
-    double mflops;
-    double mflops_best = 0;
-    double mflops_worst= 0;
-    std::vector<double> mflops_all;
-
-    ///////////////////////////////////////////////////////
-    // Set/Get the layout & grid size
-    ///////////////////////////////////////////////////////
-    int threads = GridThread::GetThreads();
-    Coordinate mpi = GridDefaultMpi(); GRID_ASSERT(mpi.size()==4);
-    Coordinate local({L,L,L,L});
-    Coordinate latt4({local[0]*mpi[0],local[1]*mpi[1],local[2]*mpi[2],local[3]*mpi[3]});
-    
-    GridCartesian         * TmpGrid   = SpaceTimeGrid::makeFourDimGrid(latt4,
-								       GridDefaultSimd(Nd,vComplex::Nsimd()),
-								       GridDefaultMpi());
-    uint64_t NP = TmpGrid->RankCount();
-    uint64_t NN = TmpGrid->NodeCount();
-    NN_global=NN;
-    uint64_t SHM=NP/NN;
-
-
-    ///////// Welcome message ////////////
-    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-    std::cout<<GridLogMessage << "Benchmark NaiveStaggered on "<<L<<"^4 local volume "<<std::endl;
-    std::cout<<GridLogMessage << "* Global volume  : "<<GridCmdVectorIntToString(latt4)<<std::endl;
-    std::cout<<GridLogMessage << "* ranks          : "<<NP  <<std::endl;
-    std::cout<<GridLogMessage << "* nodes          : "<<NN  <<std::endl;
-    std::cout<<GridLogMessage << "* ranks/node     : "<<SHM <<std::endl;
-    std::cout<<GridLogMessage << "* ranks geom     : "<<GridCmdVectorIntToString(mpi)<<std::endl;
-    std::cout<<GridLogMessage << "* Using "<<threads<<" threads"<<std::endl;
-    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-
-    ///////// Lattice Init ////////////
-    GridCartesian         * FGrid   = SpaceTimeGrid::makeFourDimGrid(latt4, GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
-    GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(FGrid);
-    
-    ///////// RNG Init ////////////
-    std::vector<int> seeds4({1,2,3,4});
-    GridParallelRNG          RNG4(FGrid);  RNG4.SeedFixedIntegers(seeds4);
-    std::cout << GridLogMessage << "Initialised RNGs" << std::endl;
-
-    RealD mass=0.1;
-    RealD c1=9.0/8.0;
-    RealD c2=-1.0/24.0;
-    RealD u0=1.0;
-
-    typedef NaiveStaggeredFermionF Action;
-    typedef typename Action::FermionField Fermion; 
-    typedef LatticeGaugeFieldF Gauge;
-    
-    Gauge Umu(FGrid);  SU<Nc>::HotConfiguration(RNG4,Umu); 
-
-    typename Action::ImplParams params;
-    Action Ds(Umu,*FGrid,*FrbGrid,mass,c1,u0,params);
-
-    ///////// Source preparation ////////////
-    Fermion src   (FGrid); random(RNG4,src);
-    Fermion src_e (FrbGrid);
-    Fermion src_o (FrbGrid);
-    Fermion r_e   (FrbGrid);
-    Fermion r_o   (FrbGrid);
-    Fermion r_eo  (FGrid);
-  
-    {
-
-      pickCheckerboard(Even,src_e,src);
-      pickCheckerboard(Odd,src_o,src);
-    
-      const int num_cases = 2;
-      std::string fmt("G/S/C ; G/O/C ; G/S/S ; G/O/S ");
-      
-      controls Cases [] = {
-	{  StaggeredKernelsStatic::OptGeneric   ,  StaggeredKernelsStatic::CommsAndCompute  ,CartesianCommunicator::CommunicatorPolicyConcurrent  },
-	{  StaggeredKernelsStatic::OptHandUnroll,  StaggeredKernelsStatic::CommsAndCompute  ,CartesianCommunicator::CommunicatorPolicyConcurrent  },
-	{  StaggeredKernelsStatic::OptInlineAsm ,  StaggeredKernelsStatic::CommsAndCompute  ,CartesianCommunicator::CommunicatorPolicyConcurrent  }
-      }; 
-
-      for(int c=0;c<num_cases;c++) {
-	
-	StaggeredKernelsStatic::Comms = Cases[c].CommsOverlap;
-	StaggeredKernelsStatic::Opt   = Cases[c].Opt;
-	CartesianCommunicator::SetCommunicatorPolicy(Cases[c].CommsAsynch);
-      
-	std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-	if ( StaggeredKernelsStatic::Opt == StaggeredKernelsStatic::OptGeneric   ) std::cout << GridLogMessage<< "* Using GENERIC Nc StaggeredKernels" <<std::endl;
-	std::cout << GridLogMessage<< "* SINGLE precision "<<std::endl;
-	std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-	
-	int nwarm = 10;
-	double t0=usecond();
-	FGrid->Barrier();
-	for(int i=0;i<nwarm;i++){
-	  Ds.DhopEO(src_o,r_e,DaggerNo);
-	}
-	FGrid->Barrier();
-	double t1=usecond();
-
-	uint64_t no    = 50;
-	uint64_t ni    = 100;
-
-	//	std::cout << GridLogMessage << " Estimate " << ncall << " calls per second"<<std::endl;
-
-	time_statistics timestat;
-	std::vector<double> t_time(no);
-	for(uint64_t i=0;i<no;i++){
-	  t0=usecond();
-	  for(uint64_t j=0;j<ni;j++){
-	    Ds.DhopEO(src_o,r_e,DaggerNo);
-	  }
-	  t1=usecond();
-	  t_time[i] = t1-t0;
-	}
-	FGrid->Barrier();
-	
-	double volume=1;  for(int mu=0;mu<Nd;mu++) volume=volume*latt4[mu];
-	double flops=((8*(3*(6+8+8)) + 7*3*2)*1.0*volume)/2;
-	double mf_hi, mf_lo, mf_err;
-	
-	timestat.statistics(t_time);
-	mf_hi = flops/timestat.min*ni;
-	mf_lo = flops/timestat.max*ni;
-	mf_err= flops/timestat.min * timestat.err/timestat.mean;
-
-	mflops = flops/timestat.mean*ni;
-	mflops_all.push_back(mflops);
-	if ( mflops_best == 0   ) mflops_best = mflops;
-	if ( mflops_worst== 0   ) mflops_worst= mflops;
-	if ( mflops>mflops_best ) mflops_best = mflops;
-	if ( mflops<mflops_worst) mflops_worst= mflops;
-	
-	std::cout<<GridLogMessage << std::fixed << std::setprecision(1)<<"Deo mflop/s =   "<< mflops << " ("<<mf_err<<") " << mf_lo<<"-"<<mf_hi <<std::endl;
-	std::cout<<GridLogMessage << std::fixed << std::setprecision(1)<<"Deo mflop/s per rank   "<< mflops/NP<<std::endl;
-	std::cout<<GridLogMessage << std::fixed << std::setprecision(1)<<"Deo mflop/s per node   "<< mflops/NN<<std::endl;
-	std::cout<<GridLogMessage << std::fixed << std::setprecision(1)<<"Deo us per call   "<< timestat.mean/ni<<std::endl;
-      
-      }
-
-      std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-      std::cout<<GridLogMessage << L<<"^4  Deo Best  mflop/s        =   "<< mflops_best << " ; " << mflops_best/NN<<" per node " <<std::endl;
-      std::cout<<GridLogMessage << L<<"^4  Deo Worst mflop/s        =   "<< mflops_worst<< " ; " << mflops_worst/NN<<" per node " <<std::endl;
-      std::cout<<GridLogMessage <<fmt << std::endl;
-      std::cout<<GridLogMessage ;
-
-      for(int i=0;i<mflops_all.size();i++){
-	std::cout<<mflops_all[i]/NN<<" ; " ;
-      }
-      std::cout<<std::endl;
-    }
-    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-    return mflops_best;
-  }
-  
  static double Clover(int L)
  {
    double mflops;
@@ -1042,7 +887,6 @@ int main (int argc, char ** argv)
  std::vector<double> clover;
  std::vector<double> dwf4;
  std::vector<double> staggered;
-  std::vector<double> naive_staggered;

  int Ls=1;
  if (do_dslash){
@@ -1070,21 +914,13 @@ int main (int argc, char ** argv)
    staggered.push_back(result);
  }

-  std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-  std::cout<<GridLogMessage << " Naive Staggered dslash 4D vectorised" <<std::endl;
-  std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-  for(int l=0;l<L_list.size();l++){
-    double result = Benchmark::NaiveStaggered(L_list[l]) ;
-    naive_staggered.push_back(result);
-  }
-

  std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
  std::cout<<GridLogMessage << " Summary table Ls="<<Ls <<std::endl;
  std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-  std::cout<<GridLogMessage << "L \t\t Clover \t\t DWF4 \t\t Staggered \t\t Naive Staggered" <<std::endl;
+  std::cout<<GridLogMessage << "L \t\t Clover \t\t DWF4 \t\t Staggered" <<std::endl;
  for(int l=0;l<L_list.size();l++){
-    std::cout<<GridLogMessage << L_list[l] <<" \t\t "<< clover[l]<<" \t\t "<<dwf4[l] << " \t\t "<< staggered[l]<<" \t\t "<<naive_staggered[l]<<std::endl;
+    std::cout<<GridLogMessage << L_list[l] <<" \t\t "<< clover[l]<<" \t\t "<<dwf4[l] << " \t\t "<< staggered[l]<<std::endl;
  }
  std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
  }
@@ -1094,14 +930,14 @@ int main (int argc, char ** argv)
    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
    std::cout<<GridLogMessage << " Per Node Summary table Ls="<<Ls <<std::endl;
    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-    std::cout<<GridLogMessage << " L \t\t Clover\t\t DWF4\t\t Staggered \t\t NaiveStag \t|\t (GF/s per node)" <<std::endl;
+    std::cout<<GridLogMessage << " L \t\t Clover\t\t DWF4\t\t Staggered (GF/s per node)" <<std::endl;
    fprintf(FP,"Per node summary table\n");
    fprintf(FP,"\n");
-    fprintf(FP,"L , Wilson, DWF4, Staggered, NaiveStag\n");
+    fprintf(FP,"L , Wilson, DWF4, Staggered, GF/s per node\n");
    fprintf(FP,"\n");
    for(int l=0;l<L_list.size();l++){
-      std::cout<<GridLogMessage << L_list[l] <<" \t\t "<< clover[l]/NN<<" \t "<<dwf4[l]/NN<< " \t "<<staggered[l]/NN<<" \t " <<naive_staggered[l]/NN<<std::endl;
-      fprintf(FP,"%d , %.0f, %.0f, %.0f, %.0f\n",L_list[l],clover[l]/NN/1000.,dwf4[l]/NN/1000.,staggered[l]/NN/1000.,naive_staggered[l]/NN/1000.);
+      std::cout<<GridLogMessage << L_list[l] <<" \t\t "<< clover[l]/NN<<" \t "<<dwf4[l]/NN<< " \t "<<staggered[l]/NN<<std::endl;
+      fprintf(FP,"%d , %.0f, %.0f, %.0f\n",L_list[l],clover[l]/NN/1000.,dwf4[l]/NN/1000.,staggered[l]/NN/1000.);
    }
    fprintf(FP,"\n");
    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
@@ -8,7 +8,6 @@ LIME=/p/home/jusers/boyle2/juwels/gm2dwf/boyle/
    --disable-gparity \
    --disable-fermion-reps \
    --with-lime=$LIME \
-    --enable-accelerator-cshift \
    --disable-unified \
    CXX=nvcc \
    LDFLAGS="-cudart shared " \
@@ -0,0 +1,104 @@
+#!/bin/bash -l
+
+# Node count scan.  Everything that could confound the comparison is held
+# fixed and only the number of Lustre clients changes:
+#
+#   local volume   8.8.32.128        = 262144 sites x 576 B = 151 MB / rank
+#   aggregation    k=2, row of 16, 8 extents of 18 MB        (identical at
+#                                                             every node count)
+#   old path       contiguous run = lLattice[0]*576 = 4.6 KB (identical too)
+#
+# So the ratio of contiguous extent between the two layout preserving paths is
+# a constant 4000x across the whole scan, and any change in their relative
+# bandwidth is a property of the client count alone.  That is the variable the
+# 4 node runs could not move.
+#
+# Trim the tail of the run list if you want a smaller allocation; the
+# geometries are independent.
+
+#SBATCH --job-name=ioScale
+#SBATCH --nodes=128
+#SBATCH --ntasks-per-node=8
+#SBATCH --cpus-per-task=7
+#SBATCH --gpus-per-node=8
+#SBATCH --time=02:00:00
+#SBATCH --account=phy157_dwf
+#SBATCH --gpu-bind=none
+#SBATCH --exclusive
+#SBATCH --mem=0
+
+cat << EOF > select_gpu
+#!/bin/bash
+export GPU_MAP=(0 1 2 3 7 6 5 4)
+export NUMA_MAP=(3 3 1 1 2 2 0 0)
+export GPU=\${GPU_MAP[\$SLURM_LOCALID]}
+export NUMA=\${NUMA_MAP[\$SLURM_LOCALID]}
+unset ROCR_VISIBLE_DEVICES
+export HIP_VISIBLE_DEVICES=\$GPU
+
+rank=\$SLURM_PROCID
+CMD="numactl -m \$NUMA -N \$NUMA  \$*"
+exec \$CMD
+EOF
+
+chmod +x ./select_gpu
+SELECT=$SLURM_SUBMIT_DIR/select_gpu
+
+root=$HOME/ParallelIO/systems/Frontier/
+source ${root}/sourceme-rocm7.2.sh
+
+export OMP_NUM_THREADS=7
+export MPICH_GPU_SUPPORT_ENABLED=1
+
+WORK=/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle/aggscale.$SLURM_JOB_ID
+mkdir -p $WORK
+cd       $WORK
+
+# The largest rung writes three files of 154 GB, so budget ~700 GB for the
+# whole scan and check the quota before submitting.
+#
+# Do NOT stripe by default.  Measured on Frontier at 128 nodes / 1024 ranks,
+# 154.6 GB record, best of 3 (MB/s):
+#
+#                     default layout    lfs setstripe -c -1 -S 8M
+#   raw                    97673                81728
+#   MPI-IO lexicographic   31675                77534
+#   aggregate             189535                86063
+#
+# The two layout preserving paths want opposite things.  Wide striping is what
+# rescues MPI-IO at scale, 32 -> 78 GB/s, which is presumably what careful
+# striping has been buying in production.  It costs the aggregate path a factor
+# of 2.2, and costs raw as well.  With the default layout MPI-IO plateaus at
+# ~30 GB/s from 64 nodes upward -- the 20-40 GB/s seen in production -- while
+# the aggregate path keeps scaling.  Uncomment only to reproduce that.
+#lfs setstripe -c -1 -S 8M $WORK
+lfs getstripe -d $WORK      # record the inherited layout; needed to read the numbers
+
+BIN=$SLURM_SUBMIT_DIR/tests/IO/Test_aggregate_io
+ARGS="--aggregate-target 4194304 --io-no-correctness --io-reps 3"
+
+run () {   # run <nodes> <grid> <mpi> <comment>
+  local nodes=$1 nranks=$(( $1 * 8 ))
+  echo
+  echo "==================================================================="
+  echo "=== SCAN nodes=$nodes ranks=$nranks  grid=$2  mpi=$3   $4"
+  echo "==================================================================="
+  srun -N$nodes -n$nranks --ntasks-per-node=8 $SELECT $BIN --grid $2 --mpi $3 \
+       --accelerator-threads 8 --shm 2048 --device-mem 32000 $ARGS
+  echo "=== exit $?"
+}
+
+#            nodes  global lattice      decomposition   record
+run    4     32.32.64.128    4.4.2.1   "4.8 GB"
+run    8     32.32.64.256    4.4.2.2   "9.7 GB"
+run   16     32.32.64.512    4.4.2.4   "19.3 GB"
+run   32     32.32.128.512   4.4.4.4   "38.6 GB"
+run   64     32.32.128.1024  4.4.4.8   "77.3 GB"
+run  128     32.32.256.1024  4.4.8.8   "154.6 GB"
+
+echo
+echo "=== done.  Files left in $WORK"
+echo "=== For cold read numbers, submit a SECOND job on a fresh allocation"
+echo "=== with WORK pointed at this directory and --io-read-only added to"
+echo "=== ARGS.  Reading back what this job wrote only measures page cache."
+ls -l $WORK
@@ -1,91 +1,76 @@
-Per node summary table
-
-L , Wilson, DWF4, Staggered, NaiveStag
-
-8 , 90, 933, 38, 23
-12 , 403, 1688, 178, 113
-16 , 188, 1647, 449, 295
-24 , 947, 1574, 674, 553
-32 , 931, 1371, 718, 643
-
 Memory Bandwidth

 Bytes, GB/s per node
-786432, 40.271620
-12582912, 433.611792
-63700992, 905.374321
-201326592, 1114.979152
-491520000, 1180.241898
+6291456, 379.297050
+100663296, 3754.674992
+509607936, 6521.472413
+1610612736, 8513.456479
+3932160000, 9018.901766
+
+
+GEMM
+
+ M, N, K, BATCH, GF/s per rank
+16, 8, 16, 256, 0.564958
+16, 16, 16, 256, 243.148058
+16, 32, 16, 256, 440.346877
+32, 8, 32, 256, 439.194136
+32, 16, 32, 256, 847.334141
+32, 32, 32, 256, 1430.892623
+64, 8, 64, 256, 1242.756741
+64, 16, 64, 256, 2196.689493
+64, 32, 64, 256, 3697.458072
+16, 8, 256, 256, 899.582627
+16, 16, 256, 256, 1673.537756
+16, 32, 256, 256, 2959.597089
+32, 8, 256, 256, 1558.858630
+32, 16, 256, 256, 2864.839445
+32, 32, 256, 256, 4810.671254
+64, 8, 256, 256, 2386.092942
+64, 16, 256, 256, 4451.665937
+64, 32, 256, 256, 5942.124095
+8, 256, 16, 256, 799.867271
+16, 256, 16, 256, 1584.624888
+32, 256, 16, 256, 1949.422338
+8, 256, 32, 256, 1389.417474
+16, 256, 32, 256, 2668.344493
+32, 256, 32, 256, 3234.162120
+8, 256, 64, 256, 2150.925128
+16, 256, 64, 256, 4012.488132
+32, 256, 64, 256, 5154.785521
+


 Communications

 Packet bytes, direction, GB/s per node
+4718592, 1, 245.026198
+4718592, 2, 251.180996
+4718592, 3, 361.110977
+4718592, 5, 247.898447
+4718592, 6, 249.867523
+4718592, 7, 359.033061
+15925248, 1, 255.030946
+15925248, 2, 264.453890
+15925248, 3, 392.949183
+15925248, 5, 256.040644
+15925248, 6, 264.681896
+15925248, 7, 392.102622
+37748736, 1, 258.823333
+37748736, 2, 268.181577
+37748736, 3, 401.478191
+37748736, 5, 258.995363
+37748736, 6, 268.206586
+37748736, 7, 400.397611


-GEMM
-
- M, N, K, BATCH, GF/s per rank fp64
-16, 8, 16, 4096, 693.316363
-16, 12, 16, 4096, 657.277058
-16, 16, 16, 4096, 711.992616
-32, 8, 32, 4096, 821.084324
-32, 12, 32, 4096, 1279.852719
-32, 16, 32, 4096, 2647.096674
-64, 8, 64, 4096, 2630.192325
-64, 12, 64, 4096, 3338.071321
-64, 16, 64, 4096, 3950.899281
-16, 8, 256, 4096, 1638.362501
-16, 12, 256, 4096, 2377.502234
-16, 16, 256, 4096, 3048.328833
-32, 8, 256, 4096, 2917.384276
-32, 12, 256, 4096, 4103.085151
-32, 16, 256, 4096, 5102.971860
-64, 8, 256, 4096, 3222.258206
-64, 12, 256, 4096, 4619.456391
-64, 16, 256, 4096, 5847.916650
-8, 256, 16, 4096, 1728.073337
-12, 256, 16, 4096, 2356.653970
-16, 256, 16, 4096, 2676.876038
-8, 256, 32, 4096, 2611.531990
-12, 256, 32, 4096, 3451.573106
-16, 256, 32, 4096, 3966.915301
-8, 256, 64, 4096, 3436.248737
-12, 256, 64, 4096, 4539.497945
-16, 256, 64, 4096, 5307.992323
-
-
-
-GEMM
-
- M, N, K, BATCH, GF/s per rank fp32
-16, 8, 16, 4096, 499.017445
-16, 12, 16, 4096, 731.543385
-16, 16, 16, 4096, 958.800786
-32, 8, 32, 4096, 1549.813550
-32, 12, 32, 4096, 2147.907502
-32, 16, 32, 4096, 2601.698596
-64, 8, 64, 4096, 3785.446233
-64, 12, 64, 4096, 5116.694843
-64, 16, 64, 4096, 6109.345016
-16, 8, 256, 4096, 1206.627737
-16, 12, 256, 4096, 1809.699599
-16, 16, 256, 4096, 2412.014053
-32, 8, 256, 4096, 2406.114488
-32, 12, 256, 4096, 3605.531907
-32, 16, 256, 4096, 4798.444037
-64, 8, 256, 4096, 4688.711196
-64, 12, 256, 4096, 6990.696301
-64, 16, 256, 4096, 9214.749925
-8, 256, 16, 4096, 2596.307289
-12, 256, 16, 4096, 3439.892562
-16, 256, 16, 4096, 3907.201036
-8, 256, 32, 4096, 3012.752067
-12, 256, 32, 4096, 3904.217583
-16, 256, 32, 4096, 4599.047092
-8, 256, 64, 4096, 3721.999042
-12, 256, 64, 4096, 5098.573927
-16, 256, 64, 4096, 6159.080872
+Per node summary table

+L , Wilson, DWF4, Staggered, GF/s per node

+8 , 155, 1386, 50
+12 , 694, 4208, 230
+16 , 1841, 6675, 609
+24 , 3934, 8573, 1641
+32 , 5083, 9771, 3086

@@ -1,3 +1,4 @@
+CLIME=`spack find --paths c-lime@2-3-9 | grep c-lime| cut -c 15-`
 ../../configure --enable-comms=mpi-auto \
 --with-lime=$CLIME \
 --enable-unified=no \
@@ -8,9 +9,8 @@
 --disable-gparity \
 --disable-fermion-reps \
 --enable-simd=GPU \
--with-gmp=$GMP \
--with-mpfr=$MPFR \
--with-openssl=$OPENSSL \
+--with-gmp=$OLCF_GMP_ROOT \
+--with-mpfr=/opt/cray/pe/gcc/mpfr/3.1.4/ \
 --disable-fermion-reps \
 CXX=hipcc MPICXX=mpicxx \
 CXXFLAGS="-fPIC -I${ROCM_PATH}/include/ -I${MPICH_DIR}/include " \
@@ -0,0 +1,14 @@
+
+echo spack
+. /autofs/nccs-svm1_home1/paboyle/spack/share/spack/setup-env.sh
+
+export CLIME=`spack find --paths c-lime | grep ^c-lime | awk '{print $2}' `
+export MPFR=`spack find --paths mpfr    | grep ^mpfr  | awk '{print $2}' `
+export OPENSSL=`spack find --paths openssl | grep openssl | awk  '{print $2}' `
+export GMP=`spack find --paths gmp      | grep ^gmp | awk '{print $2}' `
+
+module load cce/20.0.0
+module load cpe/25.09
+module load rocm/6.4.0
+export LD_LIBRARY_PATH=$CRAY_LD_LIBRARY_PATH:$LD_LIBRARY_PATH
+export LD_LIBRARY_PATH=/opt/rocm-6.4.0/lib/llvm/lib/:$LD_LIBRARY_PATH
@@ -0,0 +1,10 @@
+
+echo spack
+. /autofs/nccs-svm1_home1/paboyle/spack/share/spack/setup-env.sh
+
+
+module load cce/21.0.0
+module load cpe/26.03
+module load rocm/7.0.2
+export LD_LIBRARY_PATH=$CRAY_LD_LIBRARY_PATH:$LD_LIBRARY_PATH
+export LD_LIBRARY_PATH=/opt/rocm-7.0.2/lib/llvm/lib/:$LD_LIBRARY_PATH
@@ -0,0 +1,14 @@
+
+echo spack
+. /autofs/nccs-svm1_home1/paboyle/spack/share/spack/setup-env.sh
+
+export CLIME=`spack find --paths c-lime | grep ^c-lime | awk '{print $2}' `
+export MPFR=`spack find --paths mpfr    | grep ^mpfr  | awk '{print $2}' `
+export OPENSSL=`spack find --paths openssl | grep openssl | awk  '{print $2}' `
+export GMP=`spack find --paths gmp      | grep ^gmp | awk '{print $2}' `
+
+module load cce/21.0.0
+module load cpe/26.03
+module load rocm/7.2.0
+export LD_LIBRARY_PATH=$CRAY_LD_LIBRARY_PATH:$LD_LIBRARY_PATH
+export LD_LIBRARY_PATH=/opt/rocm-7.2.0/lib/llvm/lib/:$LD_LIBRARY_PATH
@@ -2,10 +2,6 @@
 echo spack
 . /autofs/nccs-svm1_home1/paboyle/spack/share/spack/setup-env.sh

-export CLIME=`spack find --paths c-lime | grep ^c-lime | awk '{print $2}' `
-export MPFR=`spack find --paths mpfr    | grep ^mpfr  | awk '{print $2}' `
-export OPENSSL=`spack find --paths openssl | grep openssl | awk  '{print $2}' `
-export GMP=`spack find --paths gmp      | grep ^gmp | awk '{print $2}' `

 module load cce/21.0.0
 module load cpe/26.03
@@ -1,20 +0,0 @@
-#!/usr/bin/env bash
-# Configure Grid from a build directory two levels below the source root,
-# e.g.: mkdir -p systems/mac-arm/build && cd systems/mac-arm/build && bash ../config-command-homebrew
-#
-# Prerequisites: source sourceme-homebrew.sh first.
-
-CXX=mpicxx ../../configure \
-    --enable-simd=GEN \
-    --enable-comms=mpi-auto \
-    --enable-unified=yes \
-    --prefix="${HOME}/Grid-install" \
-    --with-lime="${CLIME}" \
-    --with-openssl="${OPENSSL}" \
-    --with-gmp="${GMP}" \
-    --with-mpfr="${MPFR}" \
-    --with-fftw="${FFTW}" \
-    --enable-Sp=no \
-    --disable-fermion-reps \
-    --disable-gparity \
-    --disable-debug
@@ -1,15 +0,0 @@
-#!/usr/bin/env bash
-# Environment for building Grid on Apple Silicon Mac with Homebrew dependencies.
-# Usage: source systems/mac-arm/sourceme-homebrew.sh
-
-HOMEBREW=/opt/homebrew
-
-export GMP=${HOMEBREW}/opt/gmp
-export MPFR=${HOMEBREW}/opt/mpfr
-export FFTW=${HOMEBREW}/opt/fftw
-export OPENSSL=${HOMEBREW}/opt/openssl@3
-export CLIME=/usr/local
-
-export PATH="${HOMEBREW}/opt/open-mpi/bin:${HOMEBREW}/bin:${PATH}"
-export LDFLAGS="-L${OPENSSL}/lib"
-export CPPFLAGS="-I${OPENSSL}/include"
@@ -0,0 +1,458 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/IO/Test_aggregate_io.cc
+
+    Copyright (C) 2015
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+
+    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.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+    *************************************************************************************/
+    /*  END LEGAL */
+
+// Correctness and performance test for BINARYIO_AGGREGATE.
+//
+// Correctness, per aggregateTargetBytes:
+//   1. write via the MPI-IO lexicographic path      -> ref.bin
+//   2. write via the aggregate path                 -> agg.bin
+//   3. the two files must be byte identical         <- proves the layout matches
+//   4. read agg.bin back through the aggregate path <- proves the mirror inverts
+//   5. write and read the non-lexicographic path    -> raw.bin.  Its layout is
+//      different by construction (each rank owns one contiguous segment in
+//      rank order) so it cannot be compared byte for byte, but the NERSC and
+//      SciDAC checksums are computed from the global site index and are
+//      therefore layout independent: they must match the other two paths.
+//   6. a record written over a longer pre-existing file must leave the file at
+//      exactly offset+payload, with no trailing fragment of the old contents
+//
+// Performance: three paths, both directions, timed with the client page cache
+// dropped before every read so that a read back reports filesystem bandwidth
+// rather than memory bandwidth.  The non-lexicographic path is the zero
+// overhead reference: no transposition, no layout independence, one disjoint
+// contiguous segment per rank, which is the arrangement that reaches full
+// filesystem bandwidth on a leadership machine.  It is the upper bound the
+// other two are trying to approach.
+//
+// Options:
+//   --aggregate-target <bytes>   sweep this one target only (default: sweep
+//                                1, 1024, 64K, 4M)
+//   --io-reps <n>                repetitions in the performance section
+//                                (default 3; 0 disables it)
+//   --io-no-correctness          skip the correctness section, which reads the
+//                                whole file on one rank and is not affordable
+//                                at very large volume
+//   --io-read-only               time reads only, of files left in place by an
+//                                earlier job.  Reading back what this job just
+//                                wrote measures the client page cache; a fresh
+//                                allocation pointed at the same directory is
+//                                the only way to get a cold read without root.
+//
+// The exchange is only meaningfully exercised when the fast dimensions are
+// split across ranks; --mpi 1.1.X.Y leaves the rows of size one and the test
+// then passes vacuously.  Non-uniform AllToAllV counts additionally need an
+// odd process factor in a fast dimension and a small local volume.
+
+#include <Grid/Grid.h>
+#include <fcntl.h>
+
+using namespace Grid;
+
+/////////////////////////////////////////////////////////////////////////////
+// Compare in chunks.  Slurping both files into memory is fine for a few MB
+// and fatal for the multi-GB records this test is meant to reach.
+/////////////////////////////////////////////////////////////////////////////
+static bool FilesIdentical(std::string a,std::string b)
+{
+  std::ifstream fa(a,std::ios::binary), fb(b,std::ios::binary);
+  if ( !fa.good() || !fb.good() ) {
+    std::cout<<GridLogMessage<<"  could not open "<<a<<" and/or "<<b<<std::endl;
+    return false;
+  }
+  fa.seekg(0,std::ios::end); fb.seekg(0,std::ios::end);
+  uint64_t sa = (uint64_t)fa.tellg(), sb = (uint64_t)fb.tellg();
+  if ( sa != sb ) {
+    std::cout<<GridLogMessage<<"  size mismatch "<<sa<<" vs "<<sb<<std::endl;
+    return false;
+  }
+  fa.seekg(0,std::ios::beg); fb.seekg(0,std::ios::beg);
+
+  const uint64_t chunk = 8*1024*1024;
+  std::vector<char> va(chunk), vb(chunk);
+  uint64_t done=0;
+  while ( done < sa ) {
+    uint64_t n = std::min(chunk,sa-done);
+    fa.read(&va[0],n);
+    fb.read(&vb[0],n);
+    for(uint64_t i=0;i<n;i++){
+      if ( va[i]!=vb[i] ) {
+        std::cout<<GridLogMessage<<"  first differing byte at "<<done+i<<" of "<<sa<<std::endl;
+        return false;
+      }
+    }
+    done += n;
+  }
+  return true;
+}
+
+/////////////////////////////////////////////////////////////////////////////
+// Reading back a file we have just written measures the client page cache,
+// not the filesystem: the earlier runs of this test reported 8 GB/s on reads
+// and ~1 GB/s on writes for the same data.  POSIX_FADV_DONTNEED asks the
+// kernel to drop the cached pages for the file.  It is advisory and every
+// rank must do it, since each client caches independently, so treat this as
+// best effort rather than a guarantee of a cold read.
+/////////////////////////////////////////////////////////////////////////////
+static void DropCache(GridBase *grid,std::string file)
+{
+  grid->Barrier();
+  int fd = ::open(file.c_str(),O_RDONLY);
+  if ( fd >= 0 ) {
+#ifdef POSIX_FADV_DONTNEED
+    ::posix_fadvise(fd,0,0,POSIX_FADV_DONTNEED);
+#endif
+    ::close(fd);
+  }
+  grid->Barrier();
+}
+
+static uint64_t OptionU64(int argc,char **argv,const char *opt,uint64_t def)
+{
+  if ( GridCmdOptionExists(argv,argv+argc,opt) ) {
+    std::string arg = GridCmdOptionPayload(argv,argv+argc,opt);
+    return (uint64_t)std::stoull(arg);
+  }
+  return def;
+}
+
+int main(int argc,char **argv)
+{
+  Grid_init(&argc,&argv);
+
+  Coordinate latt = GridDefaultLatt();
+  Coordinate simd = GridDefaultSimd(Nd,vComplexD::Nsimd());
+  Coordinate mpi  = GridDefaultMpi();
+  GridCartesian grid(latt,simd,mpi);
+
+  typedef vLorentzColourMatrixD vobj;
+  typedef LorentzColourMatrixD  sobj;
+
+  GridParallelRNG pRNG(&grid);
+  pRNG.SeedFixedIntegers(std::vector<int>({1,2,3,4}));
+  LatticeGaugeFieldD Umu(&grid);
+  random(pRNG,Umu);
+
+  BinarySimpleMunger<sobj,sobj> munge;
+  const std::string format("IEEE64BIG");
+
+  const int lex = BinaryIO::BINARYIO_LEXICOGRAPHIC;
+  const int agg = BinaryIO::BINARYIO_LEXICOGRAPHIC|BinaryIO::BINARYIO_AGGREGATE;
+  const int raw = 0;   // no BINARYIO_LEXICOGRAPHIC: contiguous segment per rank
+
+  uint64_t payload = (uint64_t)grid._gsites*sizeof(sobj);
+
+  std::vector<uint64_t> targets = {1, 1024, 64*1024, 4*1024*1024};
+  if ( GridCmdOptionExists(argv,argv+argc,"--aggregate-target") ) {
+    targets.clear();
+    targets.push_back(OptionU64(argc,argv,"--aggregate-target",4*1024*1024));
+  }
+  uint64_t reps        = OptionU64(argc,argv,"--io-reps",3);
+  bool     correctness = !GridCmdOptionExists(argv,argv+argc,"--io-no-correctness");
+  // Read only: time reads of files left by an earlier job.  The only way to
+  // get a cold client cache without root is to read on an allocation that did
+  // not write the data, so run one job to write and a second, pointed at the
+  // same directory, with this flag.
+  bool     readonly    = GridCmdOptionExists(argv,argv+argc,"--io-read-only");
+  if ( readonly ) correctness = false;
+
+  std::cout<<GridLogMessage<<"Record payload "<<payload<<" bytes = "
+           <<payload/1024./1024.<<" MB, "
+           <<payload/(RealD)grid.ProcessorCount()/1024./1024.<<" MB/rank"<<std::endl;
+
+  int failures=0;
+
+  //////////////////////////////////////////////////////////////////////////
+  // Correctness
+  //////////////////////////////////////////////////////////////////////////
+  if ( correctness ) for(auto target : targets){
+
+    std::cout<<GridLogMessage<<"=== correctness, aggregateTargetBytes = "<<target<<" ==="<<std::endl;
+
+    uint32_t n1,a1,b1, n2,a2,b2, n3,a3,b3;
+    uint64_t off;
+
+    // Start from a clean slate.  The aggregate path sets the file length to
+    // exactly offset+payload; the MPI-IO path (MPI_MODE_CREATE) leaves any
+    // pre-existing tail in place.  Comparing stale files would therefore
+    // report a size mismatch that says nothing about the payload.
+    if ( grid.IsBoss() ) { ::unlink("ref.bin"); ::unlink("agg.bin"); ::unlink("raw.bin"); }
+    grid.Barrier();
+
+    off=0;
+    BinaryIO::writeLatticeObject<vobj,sobj>(Umu,"ref.bin",munge,off,format,n1,a1,b1,lex);
+
+    BinaryIO::aggregateTargetBytes = target;
+    off=0;
+    BinaryIO::writeLatticeObject<vobj,sobj>(Umu,"agg.bin",munge,off,format,n2,a2,b2,agg);
+
+    grid.Barrier();
+
+    if ( grid.IsBoss() ) {
+      if ( !FilesIdentical("ref.bin","agg.bin") ) {
+        std::cout<<GridLogError<<"  FAIL: aggregate file differs from lexicographic file"<<std::endl;
+        failures++;
+      } else {
+        std::cout<<GridLogMessage<<"  files byte identical"<<std::endl;
+      }
+    }
+    if ( (n1!=n2)||(a1!=a2)||(b1!=b2) ) {
+      std::cout<<GridLogError<<"  FAIL: checksum mismatch between paths"<<std::endl;
+      failures++;
+    }
+    // writeLatticeObject takes offset by value, so the out-parameter that
+    // IOobject sets never reaches us here and cannot be checked directly.
+    // The observable equivalent is the file length: both paths must leave the
+    // record ending at exactly offset+payload.
+    if ( grid.IsBoss() ) {
+      for(auto f : {std::string("ref.bin"),std::string("agg.bin")}){
+        std::ifstream fs(f,std::ios::binary|std::ios::ate);
+        uint64_t sz = (uint64_t)fs.tellg();
+        if ( sz != payload ) {
+          std::cout<<GridLogError<<"  FAIL: "<<f<<" is "<<sz<<" bytes, expected "<<payload<<std::endl;
+          failures++;
+        }
+      }
+    }
+
+    LatticeGaugeFieldD Uchk(&grid);
+    DropCache(&grid,"agg.bin");
+    off=0;
+    BinaryIO::readLatticeObject<vobj,sobj>(Uchk,"agg.bin",munge,off,format,n3,a3,b3,agg);
+    if ( (n3!=n1)||(a3!=a1)||(b3!=b1) ) {
+      std::cout<<GridLogError<<"  FAIL: read back checksum mismatch"<<std::endl;
+      failures++;
+    }
+    Uchk = Uchk - Umu;
+    RealD residual = norm2(Uchk);
+    std::cout<<GridLogMessage<<"  read back residual "<<residual<<std::endl;
+    if ( residual != 0.0 ) {
+      std::cout<<GridLogError<<"  FAIL: read back does not reproduce the field"<<std::endl;
+      failures++;
+    }
+
+    ////////////////////////////////////////////////////////////////////////
+    // Non-lexicographic.  Different file layout by construction, so compare
+    // by checksum and by round trip rather than by bytes.
+    ////////////////////////////////////////////////////////////////////////
+    uint32_t n4,a4,b4, n5,a5,b5;
+    off=0;
+    BinaryIO::writeLatticeObject<vobj,sobj>(Umu,"raw.bin",munge,off,format,n4,a4,b4,raw);
+    grid.Barrier();
+    if ( (n4!=n1)||(a4!=a1)||(b4!=b1) ) {
+      std::cout<<GridLogError<<"  FAIL: raw path checksum differs; the NERSC and"
+               <<" SciDAC checksums are layout independent and must agree"<<std::endl;
+      failures++;
+    }
+    DropCache(&grid,"raw.bin");
+    off=0;
+    BinaryIO::readLatticeObject<vobj,sobj>(Uchk,"raw.bin",munge,off,format,n5,a5,b5,raw);
+    if ( (n5!=n1)||(a5!=a1)||(b5!=b1) ) {
+      std::cout<<GridLogError<<"  FAIL: raw read back checksum mismatch"<<std::endl;
+      failures++;
+    }
+    Uchk = Uchk - Umu;
+    RealD rawresidual = norm2(Uchk);
+    std::cout<<GridLogMessage<<"  raw read back residual "<<rawresidual<<std::endl;
+    if ( rawresidual != 0.0 ) {
+      std::cout<<GridLogError<<"  FAIL: raw read back does not reproduce the field"<<std::endl;
+      failures++;
+    }
+  }
+
+  //////////////////////////////////////////////////////////////////////////
+  // Truncation.  offset!=0 is the case that matters: it is what ILDG and
+  // NERSC use, and it is the branch that used to only ever grow the file.
+  // The oversized starting file is made by extending a sparse one rather
+  // than writing padding from a single rank, which does not scale.
+  //////////////////////////////////////////////////////////////////////////
+  if ( correctness ) {
+    BinaryIO::aggregateTargetBytes = 4*1024*1024;
+    for(uint64_t testOffset : {(uint64_t)0, (uint64_t)1024}){
+
+      uint64_t expect = testOffset + payload;
+
+      if ( grid.IsBoss() ) {
+        { std::ofstream create("trunc.bin",std::ios::binary|std::ios::out); create.close(); }
+        int ierr = ::truncate("trunc.bin",(off_t)(expect+65536));
+        GRID_ASSERT(ierr==0);
+      }
+      grid.Barrier();
+
+      uint32_t n,a,b;
+      uint64_t off = testOffset;
+      BinaryIO::writeLatticeObject<vobj,sobj>(Umu,"trunc.bin",munge,off,format,n,a,b,agg);
+      grid.Barrier();
+
+      if ( grid.IsBoss() ) {
+        std::ifstream f("trunc.bin",std::ios::binary|std::ios::ate);
+        uint64_t sz = (uint64_t)f.tellg();
+        f.close();
+        if ( sz != expect ) {
+          std::cout<<GridLogError<<"  FAIL: offset "<<testOffset<<" left file "<<sz
+                   <<" bytes, expected "<<expect<<std::endl;
+          failures++;
+        } else {
+          std::cout<<GridLogMessage<<"  truncation ok at offset "<<testOffset
+                   <<": file is exactly "<<sz<<" bytes"<<std::endl;
+        }
+      }
+    }
+  }
+
+  //////////////////////////////////////////////////////////////////////////
+  // Performance.  Four numbers per repetition: write and read, old path and
+  // new.  Reads are preceded by a cache drop; writes are not, so a write
+  // number is "time to hand the data to the client cache and close", the
+  // same convention for both paths.
+  //////////////////////////////////////////////////////////////////////////
+  if ( reps ) {
+
+    LatticeGaugeFieldD Uio(&grid);
+    uint32_t n,a,b;
+
+    // The NERSC and SciDAC checksums are computed from the global site index,
+    // so all three layouts must produce the same values.  This costs nothing
+    // and is the only correctness check available at a volume where the byte
+    // for byte comparison (single rank, whole file) is unaffordable.
+    uint32_t cn[6],ca[6],cb[6];
+    auto agreeing = [&](const char *what,int lo,int hi){
+      for(int i=lo+1;i<=hi;i++){
+        if ( (cn[i]!=cn[lo])||(ca[i]!=ca[lo])||(cb[i]!=cb[lo]) ) {
+          std::cout<<GridLogError<<"  FAIL: "<<what<<" checksums disagree between paths"<<std::endl;
+          return false;
+        }
+      }
+      return true;
+    };
+
+    for(auto target : targets){
+
+      BinaryIO::aggregateTargetBytes = target;
+      std::cout<<GridLogMessage<<"=== performance, aggregateTargetBytes = "<<target<<" ==="<<std::endl;
+
+      std::vector<RealD> wref,wagg,wraw,rref,ragg,rraw;
+
+      for(uint64_t rep=0;rep<reps;rep++){
+
+        uint64_t off;
+
+        // Unlink only before the first repetition.  Lustre metadata cost is
+        // per file, not per byte, so rep 0 reports "create the file and write
+        // it" and the later reps report the steady state of overwriting an
+        // existing file -- which is what a multi record file does for every
+        // record after the first, and what production actually looks like.
+        if ( (rep==0) && !readonly && grid.IsBoss() ) {
+          ::unlink("ref.bin"); ::unlink("agg.bin"); ::unlink("raw.bin");
+        }
+        grid.Barrier();
+
+        if ( !readonly ) {
+
+        off=0;
+        BinaryIO::writeLatticeObject<vobj,sobj>(Umu,"ref.bin",munge,off,format,n,a,b,lex);
+        wref.push_back(BinaryIO::lastPerf.mbytesPerSecond);
+        cn[0]=n; ca[0]=a; cb[0]=b;
+
+        off=0;
+        BinaryIO::writeLatticeObject<vobj,sobj>(Umu,"agg.bin",munge,off,format,n,a,b,agg);
+        wagg.push_back(BinaryIO::lastPerf.mbytesPerSecond);
+        cn[1]=n; ca[1]=a; cb[1]=b;
+
+        off=0;
+        BinaryIO::writeLatticeObject<vobj,sobj>(Umu,"raw.bin",munge,off,format,n,a,b,raw);
+        wraw.push_back(BinaryIO::lastPerf.mbytesPerSecond);
+        cn[2]=n; ca[2]=a; cb[2]=b;
+        if ( !agreeing("write",0,2) ) failures++;
+
+        }  // !readonly
+
+        DropCache(&grid,"ref.bin");
+        off=0;
+        BinaryIO::readLatticeObject<vobj,sobj>(Uio,"ref.bin",munge,off,format,n,a,b,lex);
+        rref.push_back(BinaryIO::lastPerf.mbytesPerSecond);
+        cn[3]=n; ca[3]=a; cb[3]=b;
+
+        DropCache(&grid,"agg.bin");
+        off=0;
+        BinaryIO::readLatticeObject<vobj,sobj>(Uio,"agg.bin",munge,off,format,n,a,b,agg);
+        ragg.push_back(BinaryIO::lastPerf.mbytesPerSecond);
+        cn[4]=n; ca[4]=a; cb[4]=b;
+
+        DropCache(&grid,"raw.bin");
+        off=0;
+        BinaryIO::readLatticeObject<vobj,sobj>(Uio,"raw.bin",munge,off,format,n,a,b,raw);
+        rraw.push_back(BinaryIO::lastPerf.mbytesPerSecond);
+        cn[5]=n; ca[5]=a; cb[5]=b;
+        if ( !agreeing("read back",readonly?3:0,5) ) failures++;
+      }
+
+      if ( grid.IsBoss() ) {
+        auto report = [&](const char *name,std::vector<RealD> &v){
+          if ( v.empty() ) return;
+          RealD best=0, sum=0;
+          for(auto x : v){ if(x>best) best=x; sum+=x; }
+          // First sample includes file creation, later ones do not; quote both
+          // rather than a mean that mixes the two.
+          std::cout<<GridLogMessage<<"  PERF target="<<target<<" "<<name
+                   <<" best "<<best<<" MB/s, mean "<<sum/v.size()
+                   <<" MB/s, first(cold create) "<<v[0]<<" MB/s, samples";
+          for(auto x : v) std::cout<<" "<<x;
+          std::cout<<std::endl;
+        };
+        report("write raw       ",wraw);   // zero overhead reference
+        report("write MPI-IO    ",wref);
+        report("write aggregate ",wagg);
+        report("read  raw       ",rraw);
+        report("read  MPI-IO    ",rref);
+        report("read  aggregate ",ragg);
+
+        // Fraction of the zero overhead reference that each layout preserving
+        // path achieves.  This is the number the whole exercise is about.
+        auto best = [](std::vector<RealD> &v){ RealD m=0; for(auto x:v) if(x>m) m=x; return m; };
+        if ( !wraw.empty() && best(wraw) > 0 ) {
+          std::cout<<GridLogMessage<<"  PERF target="<<target
+                   <<" write fraction of raw: MPI-IO "<<best(wref)/best(wraw)
+                   <<"  aggregate "<<best(wagg)/best(wraw)<<std::endl;
+        }
+        if ( !rraw.empty() && best(rraw) > 0 ) {
+          std::cout<<GridLogMessage<<"  PERF target="<<target
+                   <<" read  fraction of raw: MPI-IO "<<best(rref)/best(rraw)
+                   <<"  aggregate "<<best(ragg)/best(rraw)<<std::endl;
+        }
+      }
+    }
+  }
+
+  if ( grid.IsBoss() ) {
+    if ( failures ) std::cout<<GridLogError  <<failures<<" FAILURE(S)"<<std::endl;
+    else            std::cout<<GridLogMessage<<"ALL AGGREGATE IO TESTS PASSED"<<std::endl;
+  }
+
+  Grid_finalize();
+  return failures!=0;
+}
@@ -0,0 +1,261 @@
+/*************************************************************************************
+
+ Test_fft_memory.cc
+
+ Memory growth test for PlannedFFT on a spin-colour matrix (propagator) field.
+
+ The test creates a single PlannedFFT object (which allocates FFTW plans once),
+ then repeatedly applies FFT_all_dim to the same propagator 400 times.
+
+ If PlannedFFT is working correctly the RSS should remain flat after the first
+ iteration — no new plans, no new deviceVector allocations beyond the per-call
+ pencil buffer which is freed at the end of each FFT_dim_execute call.
+
+ Build exactly like any other Grid test, e.g.:
+   make Test_fft_memory
+ or compile manually:
+   $(CXX) $(CXXFLAGS) Test_fft_memory.cc -o Test_fft_memory $(LDFLAGS)
+
+*************************************************************************************/
+
+#include <Grid/Grid.h>
+using namespace Grid;
+
+// --------------------------------------------------------------------------
+// Helper: read RSS (resident set size) in kB from /proc/self/status.
+// Returns 0 on platforms where /proc is unavailable.
+// --------------------------------------------------------------------------
+static long getCPURSSKb()
+{
+  long rss = 0;
+  FILE *fp = fopen("/proc/self/status", "r");
+  if (!fp) return -1;
+  char line[256];
+  while (fgets(line, sizeof(line), fp)) {
+    if (strncmp(line, "VmRSS:", 6) == 0) {
+      sscanf(line + 6, "%ld", &rss);
+      break;
+    }
+  }
+  fclose(fp);
+  return rss;
+}
+
+static long getGPUUsedMb()
+{
+#if defined(GRID_CUDA)
+  size_t free_bytes  = 0;
+  size_t total_bytes = 0;
+  cudaError_t err = cudaMemGetInfo(&free_bytes, &total_bytes);
+  if (err != cudaSuccess) return -1;
+  return (long)((total_bytes - free_bytes) / (1024 * 1024));
+ 
+#elif defined(GRID_HIP)
+  size_t free_bytes  = 0;
+  size_t total_bytes = 0;
+  hipError_t err = hipMemGetInfo(&free_bytes, &total_bytes);
+  if (err != hipSuccess) return -1;
+  return (long)((total_bytes - free_bytes) / (1024 * 1024));
+ 
+#else
+  return -1;   // CPU-only build: no GPU to query
+#endif
+}
+
+// ============================================================
+//  Convenience struct — one snapshot of both sides
+// ============================================================
+struct MemSnapshot {
+  long cpu_rss_kb;   // host RSS in kB  (-1 if unavailable)
+  long gpu_used_mb;  // device used in MB (-1 if no GPU)
+};
+ 
+static MemSnapshot takeSnapshot()
+{
+  MemSnapshot s;
+  s.cpu_rss_kb  = getCPURSSKb();
+  s.gpu_used_mb = getGPUUsedMb();
+  return s;
+}
+ 
+// ============================================================
+//  Pretty-print one row of the monitoring table
+// ============================================================
+static void printRow(int iter,
+                     const MemSnapshot &now,
+                     const MemSnapshot &prev)
+{
+  long cpu_delta = (now.cpu_rss_kb  >= 0 && prev.cpu_rss_kb  >= 0)
+                 ? now.cpu_rss_kb  - prev.cpu_rss_kb  : 0;
+  long gpu_delta = (now.gpu_used_mb >= 0 && prev.gpu_used_mb >= 0)
+                 ? now.gpu_used_mb - prev.gpu_used_mb : 0;
+ 
+  // Sign prefix so deltas are unambiguous
+  auto sign = [](long v) -> const char* { return v >= 0 ? "+" : ""; };
+ 
+  std::cout << GridLogMessage
+            << std::setw(6)  << iter
+            << "  CPU: " << std::setw(10) << now.cpu_rss_kb  << " kB"
+            << "  (" << sign(cpu_delta) << std::setw(7) << cpu_delta << " kB)"
+            << "  GPU: " << std::setw(7)  << now.gpu_used_mb << " MB"
+            << "  (" << sign(gpu_delta) << std::setw(5) << gpu_delta << " MB)"
+            << "\n";
+}
+
+// ============================================================
+
+int main(int argc, char **argv)
+{
+  Grid_init(&argc, &argv);
+
+  int threads = GridThread::GetThreads();
+  std::cout << GridLogMessage
+            << "Grid is setup to use " << threads << " threads" << std::endl;
+
+  // ------------------------------------------------------------------
+  // Grid setup — use whatever lattice/mpi/simd was passed on the CLI,
+  // e.g.  --grid 8.8.8.8 --mpi 1.1.1.1
+  // ------------------------------------------------------------------
+  Coordinate latt_size   = GridDefaultLatt();
+  Coordinate simd_layout = GridDefaultSimd(Nd, vComplexD::Nsimd());
+  Coordinate mpi_layout  = GridDefaultMpi();
+
+  GridCartesian GRID(latt_size, simd_layout, mpi_layout);
+
+  int vol = 1;
+  for (int d = 0; d < (int)latt_size.size(); d++) vol *= latt_size[d];
+
+  std::cout << GridLogMessage << "Lattice : ";
+  for (int d = 0; d < Nd; d++) std::cout << latt_size[d] << " ";
+  std::cout << std::endl;
+
+  // ------------------------------------------------------------------
+  // Propagator field: SpinColourMatrix = 12x12 complex, i.e.
+  // LatticePropagatorD (= Lattice<iSpinColourMatrix<vComplexD>>).
+  // This is the standard QCD quark propagator type.
+  // ------------------------------------------------------------------
+  LatticePropagatorD prop(&GRID);
+  
+  // ------------------------------------------------------------------
+  // Fill the propagator with a momentum-space plane wave,
+  // following the pattern from Test_fft.cc.
+  // We set each spin-colour component (a,b) to exp(i * sum_mu p_mu x_mu)
+  // with a fixed momentum p = (1,2,1,2).
+  // ------------------------------------------------------------------
+  Coordinate pvec({1, 2, 1, 2});
+
+  LatticeComplexD phase(&GRID);
+  LatticeComplexD coor(&GRID);
+  ComplexD ci(0.0, 1.0);
+
+  phase = Zero();
+  for (int mu = 0; mu < Nd; mu++) {
+    RealD TwoPiL = M_PI * 2.0 / latt_size[mu];
+    LatticeCoordinate(coor, mu);
+    phase = phase + (TwoPiL * pvec[mu]) * coor;
+  }
+  phase = exp(phase * ci);   // e^{i p.x}
+
+  // Broadcast the phase into every spin-colour matrix entry
+  prop = Zero();
+  prop = prop + phase;
+
+  std::cout << GridLogMessage
+            << "Propagator norm2 = " << norm2(prop) << std::endl;
+
+  // ------------------------------------------------------------------
+  // Baseline snapshot BEFORE PlannedFFT construction
+  // ------------------------------------------------------------------
+  MemSnapshot snap_before_plan = takeSnapshot();
+  std::cout << GridLogMessage
+            << "[mem] Before PlannedFFT construction"
+            << "  CPU: " << snap_before_plan.cpu_rss_kb  << " kB"
+            << "  GPU: " << snap_before_plan.gpu_used_mb << " MB"
+            << std::endl;
+ 
+  // ------------------------------------------------------------------
+  // Create the PlannedFFT — plans are allocated here ONCE for all
+  // dimensions and stored inside the object.
+  // ------------------------------------------------------------------
+  PlannedFFT<iSpinColourMatrix<vComplexD>> plannedFFT(&GRID);
+
+  // ------------------------------------------------------------------
+  // Snapshot AFTER plan construction — this is the true baseline
+  // for the loop, because cufftPlanMany itself grabs device memory.
+  // ------------------------------------------------------------------
+  MemSnapshot snap_after_plan = takeSnapshot();
+  std::cout << GridLogMessage
+            << "[mem] After  PlannedFFT construction"
+            << "  CPU: " << snap_after_plan.cpu_rss_kb  << " kB"
+            << "  GPU: " << snap_after_plan.gpu_used_mb << " MB"
+            << "  (plan overhead:"
+            << "  CPU +" << snap_after_plan.cpu_rss_kb - snap_before_plan.cpu_rss_kb << " kB"
+            << "  GPU +" << snap_after_plan.gpu_used_mb - snap_before_plan.gpu_used_mb << " MB)"
+            << std::endl;
+
+  MemoryManager::Print();
+  // ------------------------------------------------------------------
+  // 400-iteration loop.
+  // Each iteration computes the full 4d forward FFT of `prop`.
+  // We deliberately do NOT cache the result — we always start from
+  // the same `prop` so the FFT is recomputed identically each time.
+  // The point is to watch memory, not correctness.
+  // ------------------------------------------------------------------
+  const int Niter = 40;
+  const int Niter2 = 32;
+
+  // Print header for the memory table
+  std::cout << GridLogMessage
+              << "\n"
+              << std::setw(6)  << "iter"
+              << "  CPU: " << std::setw(10) << "RSS[kB]"
+              << "  (  delta  )"
+              << "  GPU: " << std::setw(7)  << "used[MB]"
+              << "  (delta)"
+              << "\n";
+ 
+  MemSnapshot snap_prev = snap_after_plan;
+
+  for (int i = 0; i < Niter; i++) {
+    std::vector<LatticePropagatorD> G;
+
+    for (int j = 0; j < Niter2; j++) {
+      LatticePropagatorD prop_fft(&GRID);
+
+      // Full 4d forward FFT using the pre-built plans
+      plannedFFT.FFT_all_dim(prop_fft, prop, FFT::forward);
+
+      G.push_back(prop_fft);
+    }
+ 
+    // cudaMemGetInfo reflects the state *after* any pooled frees have
+    // been committed, so this is accurate without an explicit sync —
+    // FFT_dim_execute already calls accelerator_barrier() internally.
+    MemSnapshot snap_now = takeSnapshot();
+    printRow(i, snap_now, snap_prev);
+    MemoryManager::Print();
+    snap_prev = snap_now;
+  }
+
+  // ------------------------------------------------------------------
+  // Summary
+  // ------------------------------------------------------------------
+  MemSnapshot snap_final = takeSnapshot();
+ 
+  long cpu_growth = snap_final.cpu_rss_kb  - snap_after_plan.cpu_rss_kb;
+  long gpu_growth = snap_final.gpu_used_mb - snap_after_plan.gpu_used_mb;
+ 
+  std::cout << GridLogMessage
+            << "\n==== Memory summary (baseline = after plan construction) ====\n"
+            << "  CPU RSS growth over " << Niter << " FFTs : "
+            << cpu_growth << " kB"
+            << (cpu_growth == 0 ? "  OK" : "  *** GROWING ***") << "\n"
+            << "  GPU used growth over " << Niter << " FFTs : "
+            << gpu_growth << " MB"
+            << (gpu_growth == 0 ? "  OK" : "  *** GROWING ***") << "\n"
+            << "  Note: first-call watermark from pool fill is expected and benign.\n"
+            << "        A leak shows as continuous growth beyond iter ~2-3.\n";
+ 
+  Grid_finalize();
+  return 0;
+}
@@ -36,8 +36,8 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 using namespace std;
 using namespace Grid;

-//gridblasHandle_t GridBLAS::gridblasHandle;
-//int            GridBLAS::gridblasInit;
+gridblasHandle_t GridBLAS::gridblasHandle;
+int            GridBLAS::gridblasInit;

 ///////////////////////
 // Tells little dirac op to use MdagM as the .Op()
@@ -0,0 +1,76 @@
+/*
+ * Isolating the hipfft HIPFFT_PARSE_ERROR on ROCm 7 / hipFFT 1.0.20.
+ *
+ * Tests three orderings with an empty rocFFT cache to find which GPU
+ * operation before plan creation triggers the failure:
+ *   A) hipMalloc only            — hypothesis: passes (no async GPU work)
+ *   B) hipMalloc + hipMemset     — hypothesis: fails  (async GPU work in flight)
+ *   C) hipMalloc + hipMemset     — hypothesis: passes (work completed before plan)
+ *      + hipDeviceSynchronize
+ *
+ * Compile:
+ *   hipcc -o Test_hipfft_bug_fail Test_hipfft_bug_fail.cc -lhipfft
+ *
+ * Run with empty cache:
+ *   rm -rf ~/.cache/
+ *   ./Test_hipfft_bug_fail
+ */
+
+#include <cstdio>
+#include <hipfft/hipfft.h>
+#include <hip/hip_runtime.h>
+
+static const char *res(hipfftResult rv) {
+  return rv == HIPFFT_SUCCESS ? "SUCCESS" : "PARSE_ERROR";
+}
+
+static hipfftResult makePlan(int G, int howmany) {
+  int n[] = {G};
+  hipfftHandle p;
+  size_t workSize = 0;
+  hipfftCreate(&p);
+  hipfftResult rv = hipfftMakePlanMany(p, 1, n,
+                                       nullptr, 1, G, nullptr, 1, G,
+                                       HIPFFT_Z2Z, howmany, &workSize);
+  hipfftDestroy(p);
+  return rv;
+}
+
+int main(void) {
+  hipDeviceProp_t prop;
+  hipGetDeviceProperties(&prop, 0);
+  printf("Device: %s\n", prop.name);
+#ifdef hipfftVersionMinor
+  printf("hipFFT version: %d.%d.%d\n\n",
+         hipfftVersionMajor, hipfftVersionMinor, hipfftVersionPatch);
+#endif
+
+  for (int G : {4, 8, 16, 32}) {
+    int howmany = 512;
+    long nelems = (long)G * howmany;
+    hipfftDoubleComplex *buf = nullptr;
+    hipMalloc(&buf, nelems * sizeof(hipfftDoubleComplex));
+
+    // Tests ordered so each runs before a prior success can populate the cache.
+
+    // B first: hipMalloc + hipMemset (async GPU work in flight)
+    // If this fails, A (no hipMemset) will pass, confirming hipMemset is the trigger.
+    hipMemset(buf, 0, nelems * sizeof(hipfftDoubleComplex));
+    hipfftResult rvB = makePlan(G, howmany);
+    printf("G=%-4d  B) hipMalloc + hipMemset       : %s\n", G, res(rvB));
+
+    // C: hipMalloc + hipMemset + sync — does syncing before plan creation fix it?
+    hipMemset(buf, 0, nelems * sizeof(hipfftDoubleComplex));
+    hipDeviceSynchronize();
+    hipfftResult rvC = makePlan(G, howmany);
+    printf("G=%-4d  C) hipMalloc + hipMemset + sync: %s\n", G, res(rvC));
+
+    // A last: hipMalloc only, no async GPU work — should always pass
+    hipfftResult rvA = makePlan(G, howmany);
+    printf("G=%-4d  A) hipMalloc only             : %s\n\n", G, res(rvA));
+
+    hipFree(buf);
+  }
+
+  return 0;
+}
@@ -0,0 +1,61 @@
+/*
+ * Minimal program demonstrating the workaround for the hipfft ROCm 7 bug.
+ *
+ * Workaround: create the hipfft plan BEFORE any hipMalloc.  Plan creation
+ * for G < 32 then succeeds even with an empty rocFFT cache.
+ *
+ * Compile:
+ *   hipcc -o Test_hipfft_bug_pass Test_hipfft_bug_pass.cc -lhipfft
+ *
+ * Run:
+ *   rm -rf ~/.cache/rocfft
+ *   ./Test_hipfft_bug_pass
+ *
+ * Expected: all G values succeed.
+ * Compare with Test_hipfft_bug_fail.cc which uses the opposite ordering.
+ */
+
+#include <cstdio>
+#include <hipfft/hipfft.h>
+#include <hip/hip_runtime.h>
+
+int main(void) {
+  hipDeviceProp_t prop;
+  hipGetDeviceProperties(&prop, 0);
+  printf("Device: %s\n", prop.name);
+#ifdef hipfftVersionMinor
+  printf("hipFFT version: %d.%d.%d\n\n",
+         hipfftVersionMajor, hipfftVersionMinor, hipfftVersionPatch);
+#endif
+
+  for (int G : {8, 16, 32}) {
+    int howmany = 512;
+    int n[] = {G};
+    long nelems = (long)G * howmany;
+
+    // Plan created BEFORE hipMalloc — succeeds for all G
+    hipfftHandle p;
+    size_t workSize = 0;
+    hipfftCreate(&p);
+    hipfftResult rv = hipfftMakePlanMany(p, 1, n,
+                                         nullptr, 1, G, nullptr, 1, G,
+                                         HIPFFT_Z2Z, howmany, &workSize);
+    printf("G=%-4d  plan-then-hipMalloc: %d (%s)\n",
+           G, (int)rv, rv == HIPFFT_SUCCESS ? "HIPFFT_SUCCESS" : "HIPFFT_PARSE_ERROR");
+
+    if (rv == HIPFFT_SUCCESS) {
+      hipfftDoubleComplex *buf = nullptr;
+      hipMalloc(&buf, nelems * sizeof(hipfftDoubleComplex));
+      hipMemset(buf, 0, nelems * sizeof(hipfftDoubleComplex));
+      rv = hipfftExecZ2Z(p, buf, buf, HIPFFT_FORWARD);
+      hipDeviceSynchronize();
+      printf("G=%-4d  execFwd:             %d (%s)\n",
+             G, (int)rv, rv == HIPFFT_SUCCESS ? "HIPFFT_SUCCESS" : "FAILED");
+      hipFree(buf);
+    }
+
+    hipfftDestroy(p);
+  }
+
+  return 0;
+}
@@ -0,0 +1,142 @@
+/*
+ * Minimal reproducer for hipfftMakePlanMany / hipfftPlanMany failures.
+ *
+ * Compile on Frontier (no Grid headers needed):
+ *   hipcc -o Test_hipfft_minimal Test_hipfft_minimal.cc -lhipfft
+ *
+ * Run:
+ *   ./Test_hipfft_minimal
+ */
+
+#include <cstdio>
+#include <cstdlib>
+#include <hipfft/hipfft.h>
+#include <hip/hip_runtime.h>
+
+static const char *hipfftResultString(hipfftResult r) {
+  switch (r) {
+  case HIPFFT_SUCCESS:                  return "HIPFFT_SUCCESS";
+  case HIPFFT_INVALID_PLAN:             return "HIPFFT_INVALID_PLAN";
+  case HIPFFT_ALLOC_FAILED:             return "HIPFFT_ALLOC_FAILED";
+  case HIPFFT_INVALID_TYPE:             return "HIPFFT_INVALID_TYPE";
+  case HIPFFT_INVALID_VALUE:            return "HIPFFT_INVALID_VALUE";
+  case HIPFFT_INTERNAL_ERROR:           return "HIPFFT_INTERNAL_ERROR";
+  case HIPFFT_EXEC_FAILED:              return "HIPFFT_EXEC_FAILED";
+  case HIPFFT_SETUP_FAILED:             return "HIPFFT_SETUP_FAILED";
+  case HIPFFT_INVALID_SIZE:             return "HIPFFT_INVALID_SIZE";
+  case HIPFFT_UNALIGNED_DATA:           return "HIPFFT_UNALIGNED_DATA";
+  case HIPFFT_INCOMPLETE_PARAMETER_LIST:return "HIPFFT_INCOMPLETE_PARAMETER_LIST";
+  case HIPFFT_INVALID_DEVICE:           return "HIPFFT_INVALID_DEVICE";
+  case HIPFFT_PARSE_ERROR:              return "HIPFFT_PARSE_ERROR";
+  case HIPFFT_NO_WORKSPACE:             return "HIPFFT_NO_WORKSPACE";
+  case HIPFFT_NOT_IMPLEMENTED:          return "HIPFFT_NOT_IMPLEMENTED";
+  case HIPFFT_NOT_SUPPORTED:            return "HIPFFT_NOT_SUPPORTED";
+  default:                              return "UNKNOWN";
+  }
+}
+
+// Plan creation + execution for (G, howmany).
+// Tests two orderings to isolate whether a prior hipMalloc poisons hipfft
+// plan creation for small G on ROCm 7:
+//   A) plan BEFORE hipMalloc  — hypothesis: succeeds
+//   B) hipMalloc BEFORE plan  — hypothesis: fails for G < 32
+static void tryPlanAndExec(int G, long howmany) {
+  int n[] = {G};
+  long nelems = (long)G * howmany;
+
+  printf("--- G=%-4d  howmany=%-10ld  total_elems=%-12ld ---\n",
+         G, howmany, nelems);
+
+  // --- A: create plan first, allocate buffer afterwards ---
+  {
+    hipfftHandle p;
+    size_t workSize = 0;
+    hipfftCreate(&p);
+    hipfftResult rv = hipfftMakePlanMany(p, 1, n,
+                                         nullptr, 1, G, nullptr, 1, G,
+                                         HIPFFT_Z2Z, (int)howmany, &workSize);
+    printf("  plan-first  create : %d (%s)\n", (int)rv, hipfftResultString(rv));
+    if (rv == HIPFFT_SUCCESS) {
+      hipfftDoubleComplex *buf = nullptr;
+      hipMalloc(&buf, nelems * sizeof(hipfftDoubleComplex));
+      hipMemset(buf, 0, nelems * sizeof(hipfftDoubleComplex));
+      rv = hipfftExecZ2Z(p, buf, buf, HIPFFT_FORWARD);
+      hipDeviceSynchronize();
+      printf("  plan-first  execFwd: %d (%s)\n", (int)rv, hipfftResultString(rv));
+      hipFree(buf);
+    }
+    hipfftDestroy(p);
+  }
+
+  // --- B: hipMalloc first, create plan afterwards ---
+  {
+    hipfftDoubleComplex *buf = nullptr;
+    hipMalloc(&buf, nelems * sizeof(hipfftDoubleComplex));
+    hipMemset(buf, 0, nelems * sizeof(hipfftDoubleComplex));
+
+    hipfftHandle p;
+    size_t workSize = 0;
+    hipfftCreate(&p);
+    hipfftResult rv = hipfftMakePlanMany(p, 1, n,
+                                         nullptr, 1, G, nullptr, 1, G,
+                                         HIPFFT_Z2Z, (int)howmany, &workSize);
+    printf("  malloc-first create : %d (%s)\n", (int)rv, hipfftResultString(rv));
+    if (rv == HIPFFT_SUCCESS) {
+      rv = hipfftExecZ2Z(p, buf, buf, HIPFFT_FORWARD);
+      hipDeviceSynchronize();
+      printf("  malloc-first execFwd: %d (%s)\n", (int)rv, hipfftResultString(rv));
+    }
+    hipfftDestroy(p);
+    hipFree(buf);
+  }
+
+  printf("\n");
+}
+
+int main(void) {
+  // Print HIP device info
+  int device = 0;
+  hipGetDevice(&device);
+  hipDeviceProp_t prop;
+  hipGetDeviceProperties(&prop, device);
+  printf("Device %d: %s  warpSize=%d\n\n", device, prop.name, prop.warpSize);
+
+#ifdef hipfftVersionMinor
+  printf("hipFFT version: %d.%d.%d\n\n",
+         hipfftVersionMajor, hipfftVersionMinor, hipfftVersionPatch);
+#endif
+
+  // Original sweep with small howmany (these passed first time)
+  printf("=== Small howmany (original sweep) ===\n\n");
+  for (int G : {4, 8, 12, 16, 24, 32, 48, 64})
+    tryPlanAndExec(G, 512);
+
+  // Grid-realistic howmany values derived from actual lattice geometries.
+  // howmany = Ncomp * product(ldimensions[d] for d != dim)
+  // For LatticeComplexD: Ncomp=1.
+  printf("=== Grid-realistic parameters ===\n\n");
+
+  // --grid 16.16.16.16  4D FFT  (KNOWN TO FAIL in Grid)
+  // Each dim: G=16, Nperp=16^3=4096
+  tryPlanAndExec(16, 4096);
+
+  // --grid 32.32.32.32  4D FFT  (KNOWN TO SUCCEED in Grid)
+  // Each dim: G=32, Nperp=32^3=32768
+  tryPlanAndExec(32, 32768);
+
+  // --grid 32.32.32.32 Ls=8  5D DWF FFT  (KNOWN TO FAIL on dim 0 in Grid)
+  // dim 0: G=8,  Nperp=32^4=1048576
+  tryPlanAndExec(8, 1048576);
+  // dim 1-4: G=32, Nperp=8*32^3=262144
+  tryPlanAndExec(32, 262144);
+
+  // Extra intermediate cases to bracket the failure
+  tryPlanAndExec(16, 1024);
+  tryPlanAndExec(16, 2048);
+  tryPlanAndExec(16, 8192);
+  tryPlanAndExec(8,  4096);
+  tryPlanAndExec(8,  65536);
+  tryPlanAndExec(8,  262144);
+
+  return 0;
+}
@@ -0,0 +1,168 @@
+/*
+ * Reproducer for HIPFFT_PARSE_ERROR (error 12) from hipfftMakePlanMany on
+ * ROCm 7 / hipFFT 1.0.20 (Frontier, MI210 login and MI250X compute nodes).
+ *
+ * Observed failure: G < 32 returns HIPFFT_PARSE_ERROR from all three plan
+ * creation APIs (hipfftPlanMany, hipfftMakePlanMany, hipfftPlan1d) when a
+ * device buffer is allocated and zeroed with hipMalloc+hipMemset before the
+ * plan creation call.  G >= 32 succeeds.
+ *
+ * Contrast with Test_hipfft_minimal.cc (plan-first ordering) which passes
+ * for all G even with an empty rocFFT cache.
+ *
+ * Compile on Frontier (no Grid headers needed):
+ *   hipcc -o Test_hipfft_repro Test_hipfft_repro.cc -lhipfft
+ *
+ * Run with empty cache to reproduce the failure:
+ *   rm -rf ~/.cache/rocfft
+ *   ./Test_hipfft_repro
+ */
+
+#include <cstdio>
+#include <cstdlib>
+#include <hipfft/hipfft.h>
+#include <hip/hip_runtime.h>
+
+static const char *hipfftResultString(hipfftResult r) {
+  switch (r) {
+  case HIPFFT_SUCCESS:                  return "HIPFFT_SUCCESS";
+  case HIPFFT_INVALID_PLAN:             return "HIPFFT_INVALID_PLAN";
+  case HIPFFT_ALLOC_FAILED:             return "HIPFFT_ALLOC_FAILED";
+  case HIPFFT_INVALID_TYPE:             return "HIPFFT_INVALID_TYPE";
+  case HIPFFT_INVALID_VALUE:            return "HIPFFT_INVALID_VALUE";
+  case HIPFFT_INTERNAL_ERROR:           return "HIPFFT_INTERNAL_ERROR";
+  case HIPFFT_EXEC_FAILED:              return "HIPFFT_EXEC_FAILED";
+  case HIPFFT_SETUP_FAILED:             return "HIPFFT_SETUP_FAILED";
+  case HIPFFT_INVALID_SIZE:             return "HIPFFT_INVALID_SIZE";
+  case HIPFFT_UNALIGNED_DATA:           return "HIPFFT_UNALIGNED_DATA";
+  case HIPFFT_INCOMPLETE_PARAMETER_LIST:return "HIPFFT_INCOMPLETE_PARAMETER_LIST";
+  case HIPFFT_INVALID_DEVICE:           return "HIPFFT_INVALID_DEVICE";
+  case HIPFFT_PARSE_ERROR:              return "HIPFFT_PARSE_ERROR";
+  case HIPFFT_NO_WORKSPACE:             return "HIPFFT_NO_WORKSPACE";
+  case HIPFFT_NOT_IMPLEMENTED:          return "HIPFFT_NOT_IMPLEMENTED";
+  case HIPFFT_NOT_SUPPORTED:            return "HIPFFT_NOT_SUPPORTED";
+  default:                              return "UNKNOWN";
+  }
+}
+
+// Plan creation + execution for (G, howmany) using hipfftCreate+hipfftMakePlanMany.
+// This is the path Grid's FFT.h now uses.
+static void tryPlanAndExec(int G, long howmany) {
+  int n[] = {G};
+  long nelems = (long)G * howmany;
+
+  printf("--- G=%-4d  howmany=%-10ld  total_elems=%-12ld ---\n",
+         G, howmany, nelems);
+
+  // Allocate device buffer (hipfftDoubleComplex = 16 bytes each)
+  hipfftDoubleComplex *dbuf = nullptr;
+  hipError_t herr = hipMalloc(&dbuf, nelems * sizeof(hipfftDoubleComplex));
+  if (herr != hipSuccess) {
+    printf("  hipMalloc failed (%d) for %ld elems — skipping\n\n", (int)herr, nelems);
+    return;
+  }
+  hipMemset(dbuf, 0, nelems * sizeof(hipfftDoubleComplex));
+
+  // 1. hipfftPlanMany (one-step, nullptr embed) — current Grid path
+  {
+    hipfftHandle p;
+    hipfftResult rv = hipfftPlanMany(&p, 1, n,
+                                     nullptr, 1, G,
+                                     nullptr, 1, G,
+                                     HIPFFT_Z2Z, (int)howmany);
+    printf("  hipfftPlanMany   create : %d (%s)\n", (int)rv, hipfftResultString(rv));
+    if (rv == HIPFFT_SUCCESS) {
+      rv = hipfftExecZ2Z(p, dbuf, dbuf, HIPFFT_FORWARD);
+      hipDeviceSynchronize();
+      printf("  hipfftPlanMany   execFwd: %d (%s)\n", (int)rv, hipfftResultString(rv));
+      hipfftDestroy(p);
+    }
+  }
+
+  // 2. hipfftCreate + hipfftMakePlanMany (two-step) — also current Grid path
+  {
+    hipfftHandle p;
+    size_t workSize = 0;
+    hipfftResult rc = hipfftCreate(&p);
+    if (rc == HIPFFT_SUCCESS) {
+      hipfftResult rv = hipfftMakePlanMany(p, 1, n,
+                                           nullptr, 1, G,
+                                           nullptr, 1, G,
+                                           HIPFFT_Z2Z, (int)howmany, &workSize);
+      printf("  hipfftMakePlanMany      : %d (%s)  workSize=%zu\n",
+             (int)rv, hipfftResultString(rv), workSize);
+      if (rv == HIPFFT_SUCCESS) {
+        rv = hipfftExecZ2Z(p, dbuf, dbuf, HIPFFT_FORWARD);
+        hipDeviceSynchronize();
+        printf("  hipfftMakePlanMany exec : %d (%s)\n", (int)rv, hipfftResultString(rv));
+      }
+      hipfftDestroy(p);
+    } else {
+      printf("  hipfftCreate            : %d (%s)\n", (int)rc, hipfftResultString(rc));
+    }
+  }
+
+  // 3. hipfftPlan1d (simplest API, batch = howmany)
+  {
+    hipfftHandle p;
+    hipfftResult rv = hipfftPlan1d(&p, G, HIPFFT_Z2Z, (int)howmany);
+    printf("  hipfftPlan1d     create : %d (%s)\n", (int)rv, hipfftResultString(rv));
+    if (rv == HIPFFT_SUCCESS) {
+      rv = hipfftExecZ2Z(p, dbuf, dbuf, HIPFFT_FORWARD);
+      hipDeviceSynchronize();
+      printf("  hipfftPlan1d     execFwd: %d (%s)\n", (int)rv, hipfftResultString(rv));
+      hipfftDestroy(p);
+    }
+  }
+
+  hipFree(dbuf);
+  printf("\n");
+}
+
+int main(void) {
+  // Print HIP device info
+  int device = 0;
+  hipGetDevice(&device);
+  hipDeviceProp_t prop;
+  hipGetDeviceProperties(&prop, device);
+  printf("Device %d: %s  warpSize=%d\n\n", device, prop.name, prop.warpSize);
+
+#ifdef hipfftVersionMinor
+  printf("hipFFT version: %d.%d.%d\n\n",
+         hipfftVersionMajor, hipfftVersionMinor, hipfftVersionPatch);
+#endif
+
+  // Original sweep with small howmany (these passed first time)
+  printf("=== Small howmany (original sweep) ===\n\n");
+  for (int G : {4, 8, 12, 16, 24, 32, 48, 64})
+    tryPlanAndExec(G, 512);
+
+  // Grid-realistic howmany values derived from actual lattice geometries.
+  // howmany = Ncomp * product(ldimensions[d] for d != dim)
+  // For LatticeComplexD: Ncomp=1.
+  printf("=== Grid-realistic parameters ===\n\n");
+
+  // --grid 16.16.16.16  4D FFT  (KNOWN TO FAIL in Grid)
+  // Each dim: G=16, Nperp=16^3=4096
+  tryPlanAndExec(16, 4096);
+
+  // --grid 32.32.32.32  4D FFT  (KNOWN TO SUCCEED in Grid)
+  // Each dim: G=32, Nperp=32^3=32768
+  tryPlanAndExec(32, 32768);
+
+  // --grid 32.32.32.32 Ls=8  5D DWF FFT  (KNOWN TO FAIL on dim 0 in Grid)
+  // dim 0: G=8,  Nperp=32^4=1048576
+  tryPlanAndExec(8, 1048576);
+  // dim 1-4: G=32, Nperp=8*32^3=262144
+  tryPlanAndExec(32, 262144);
+
+  // Extra intermediate cases to bracket the failure
+  tryPlanAndExec(16, 1024);
+  tryPlanAndExec(16, 2048);
+  tryPlanAndExec(16, 8192);
+  tryPlanAndExec(8,  4096);
+  tryPlanAndExec(8,  65536);
+  tryPlanAndExec(8,  262144);
+
+  return 0;
+}
@@ -1,362 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid
-
-    Source file: tests/debug/Test_staggered_hdcg.cc
-
-    Authors: Thomas Blum, Peter Boyle
-
-    HDCG (Hierarchical Deflation Conjugate Gradient) multigrid solver
-    for naive staggered fermions, based on arXiv:2409.03904.
-
-    Adapts the DWF HDCG infrastructure (Test_general_coarse_hdcg_phys48.cc) to:
-      - NaiveStaggeredFermion (nearest-neighbour only, no Naik 3-hop term)
-      - 4D SchurStaggeredOperator:  Mpc = m^2 - D_oe * D_eo  (hermitian, positive-definite)
-      - vColourVector fine field type (staggered has colour but no spin)
-      - NextToNearestStencilGeometry4D: 33-point coarse stencil
-
-    Stencil count: D_oe*D_eo has 2-hop fine range.  With blocking B >= 2 the coarse
-    shifts have L1-distance <= 2, giving 33 stencil points in 4D:
-      1 (identity) + 8 (+-e_mu) + 24 (+-e_mu +- e_nu).
-    NaiveStaggeredFermion has no Naik term, so any B >= 2 suffices.
-    To extend to ImprovedStaggeredFermion later, use B >= 6.
-
-    Reference: arXiv:2409.03904 (mrhs hermitian multigrid for DWF).
-
-    Usage (after build):
-      ./Test_staggered_hdcg --grid 16.16.16.16 --mpi 1.1.1.1
-
-*************************************************************************************/
-#include <Grid/Grid.h>
-#include <Grid/algorithms/iterative/AdefMrhs.h>
-#include <Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h>
-#include <Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczosCoarse.h>
-
-using namespace Grid;
-
-// Non-converging CG used as a smoother (fixed number of iterations)
-template<class Field>
-class CGSmoother : public LinearFunction<Field>
-{
-public:
-  typedef LinearOperatorBase<Field> FineOperator;
-  FineOperator &_Op;
-  int iters;
-  CGSmoother(int _iters, FineOperator &Op) : _Op(Op), iters(_iters) {}
-  void operator()(const Field &in, Field &out)
-  {
-    ConjugateGradient<Field> CG(0.0, iters, false);
-    out = Zero();
-    CG(_Op, in, out);
-  }
-};
-
-int main(int argc, char **argv)
-{
-  fprintf(stderr, "TRACE: entering main\n"); fflush(stderr);
-  Grid_init(&argc, &argv);
-  fprintf(stderr, "TRACE: Grid_init done\n"); fflush(stderr);
-
-  //--------------------------------------------------------------------
-  // Parameters — tune for production
-  //--------------------------------------------------------------------
-  const int nbasis = 24;   // near-null space dimension
-  const int cb     = 0;    // even checkerboard
-
-  RealD mass = 0.05;
-
-  // NaiveStaggeredFermion: nearest-neighbour hop only (no Naik term).
-  // c1 = coefficient of the hopping term (1.0 = standard normalisation).
-  // u0 = tadpole factor (1.0 = no tadpole improvement).
-  RealD c1 = 1.0;
-  RealD u0 = 1.0;
-
-  //--------------------------------------------------------------------
-  // Grids
-  // Fine:   UGrid (4D full), UrbGrid (4D red-black)
-  // Coarse: Coarse4d  with dimensions = GridDefaultLatt() / Block
-  //
-  // Recommended: GridDefaultLatt() >= 16^4, Block = {4,4,4,4}
-  // NaiveStaggeredFermion works with any Block >= {2,2,2,2}
-  //--------------------------------------------------------------------
-  fprintf(stderr, "TRACE: making UGrid\n"); fflush(stderr);
-  GridCartesian *UGrid = SpaceTimeGrid::makeFourDimGrid(
-      GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi());
-  fprintf(stderr, "TRACE: making UrbGrid\n"); fflush(stderr);
-  GridRedBlackCartesian *UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
-
-  Coordinate Block({4, 4, 4, 4});
-  Coordinate clatt = GridDefaultLatt();
-  for (int d = 0; d < clatt.size(); d++) clatt[d] /= Block[d];
-  Coordinate csimd = GridDefaultSimd(Nd, vComplex::Nsimd());
-  Coordinate cmpi  = GridDefaultMpi();
-  fprintf(stderr, "TRACE: making Coarse4d clatt=%d %d %d %d simd=%d %d %d %d mpi=%d %d %d %d Nsimd=%d\n",
-          clatt[0],clatt[1],clatt[2],clatt[3],
-          csimd[0],csimd[1],csimd[2],csimd[3],
-          cmpi[0],cmpi[1],cmpi[2],cmpi[3],
-          (int)vComplex::Nsimd()); fflush(stderr);
-
-  GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, csimd, cmpi);
-  fprintf(stderr, "TRACE: Coarse4d made\n"); fflush(stderr);
-
-  //--------------------------------------------------------------------
-  // RNG + gauge field
-  //--------------------------------------------------------------------
-  fprintf(stderr, "TRACE: RNG4\n"); fflush(stderr);
-  GridParallelRNG RNG4(UGrid);    RNG4.SeedFixedIntegers({1,2,3,4});
-  fprintf(stderr, "TRACE: RNGrb\n"); fflush(stderr);
-  GridParallelRNG RNGrb(UGrid); RNGrb.SeedFixedIntegers({5,6,7,8}); // must use full grid, not UrbGrid
-  fprintf(stderr, "TRACE: Umu\n"); fflush(stderr);
-  LatticeGaugeField Umu(UGrid);
-  fprintf(stderr, "TRACE: HotConfig\n"); fflush(stderr);
-  SU<Nc>::HotConfiguration(RNG4, Umu);
-  fprintf(stderr, "TRACE: NaiveStaggeredFermionD\n"); fflush(stderr);
-  NaiveStaggeredFermionD Ds(Umu, *UGrid, *UrbGrid, mass, c1, u0);
-  fprintf(stderr, "TRACE: SchurStaggeredOperator\n"); fflush(stderr);
-  SchurStaggeredOperator<NaiveStaggeredFermionD, LatticeStaggeredFermionD> HermOp(Ds);
-  fprintf(stderr, "TRACE: HermOp done\n"); fflush(stderr);
-
-  //--------------------------------------------------------------------
-  // Subspace: inverse-iteration near-null vectors
-  //
-  // CreateSubspace applies CG (4 solves, tol=1e-4) to random noise vectors,
-  // converging naturally to the low modes of HermOp without needing spectral
-  // bound tuning.  Switch to CreateSubspaceChebyshevNew once the spectrum is
-  // well characterised (hi ~ 5.0 for naive staggered SchurStaggeredOperator).
-  //--------------------------------------------------------------------
-  typedef Aggregation<vColourVector, vTComplex, nbasis> Subspace;
-  Subspace Aggregates(Coarse4d, UrbGrid, cb);
-
-  Aggregates.CreateSubspace(RNGrb, HermOp);
-  Aggregates.Orthogonalise();
-
-  //--------------------------------------------------------------------
-  // Coarse geometry: NextToNearestStencilGeometry4D
-  //   hops=2  ->  33 stencil points in 4D
-  //--------------------------------------------------------------------
-  NextToNearestStencilGeometry4D geom(Coarse4d);
-
-  std::cout << GridLogMessage << "Coarse stencil: " << geom.npoint << " points" << std::endl;
-
-  //--------------------------------------------------------------------
-  // Single-RHS coarse operator (used for correctness check below)
-  //--------------------------------------------------------------------
-  typedef GeneralCoarsenedMatrix<vColourVector, vTComplex, nbasis> LittleDiracOp;
-  typedef LittleDiracOp::CoarseVector CoarseVector;
-
-  LittleDiracOp LDO(geom, UrbGrid, Coarse4d);
-  LDO.CoarsenOperator(HermOp, Aggregates);
-
-  //--------------------------------------------------------------------
-  // Correctness check: P M_fine P^T c  ≈  M_coarse c
-  //
-  // Promote a random coarse vector into the fine subspace, apply the
-  // fine operator, project back, and compare with the coarse operator
-  // applied directly.  Error should be at the level of subspace
-  // approximation quality (smaller = better basis vectors).
-  //--------------------------------------------------------------------
-  {
-    GridParallelRNG RNGc(Coarse4d); RNGc.SeedFixedIntegers({9,10,11,12});
-    CoarseVector c_src(Coarse4d), c_ldop(Coarse4d), c_proj(Coarse4d);
-    random(RNGc, c_src);
-
-    LatticeStaggeredFermionD f_v(UrbGrid), f_Mv(UrbGrid);
-    Aggregates.PromoteFromSubspace(c_src, f_v);
-    HermOp.Op(f_v, f_Mv);
-    Aggregates.ProjectToSubspace(c_proj, f_Mv);
-
-    LDO.M(c_src, c_ldop);
-
-    c_proj -= c_ldop;
-    RealD err = norm2(c_proj) / norm2(c_ldop);
-    std::cout << GridLogMessage
-              << "Coarsen check |P*M_fine - M_coarse| / |M_coarse| = " << err << std::endl;
-  }
-
-  //--------------------------------------------------------------------
-  // Multi-RHS coarse grid
-  //
-  // The extra leading dimension holds nrhs right-hand sides packed into
-  // SIMD lanes, matching the pattern of Test_general_coarse_hdcg_phys48.
-  //--------------------------------------------------------------------
-  const int nrhs = vComplex::Nsimd() * 2;
-
-  Coordinate mpi   = GridDefaultMpi();
-  Coordinate rhMpi ({1,    mpi[0],  mpi[1],  mpi[2],  mpi[3]});
-  Coordinate rhLatt({nrhs, clatt[0], clatt[1], clatt[2], clatt[3]});
-  Coordinate rhSimd({vComplex::Nsimd(), 1, 1, 1, 1});
-
-  GridCartesian *CoarseMrhs = new GridCartesian(rhLatt, rhSimd, rhMpi);
-
-  typedef MultiGeneralCoarsenedMatrix<vColourVector, vTComplex, nbasis> MultiCoarseOp;
-  MultiCoarseOp mrhs(geom, CoarseMrhs);
-  mrhs.CoarsenOperator(HermOp, Aggregates, Coarse4d);
-
-  //--------------------------------------------------------------------
-  // Coarse-grid Lanczos for deflation
-  //--------------------------------------------------------------------
-  typedef HermitianLinearOperator<MultiCoarseOp, CoarseVector> MrhsHermOp;
-  MrhsHermOp MrhsCoarseOp(mrhs);
-
-  // Estimate spectral bounds for Lanczos Chebyshev filter
-  CoarseVector pm_src(CoarseMrhs); pm_src = ComplexD(1.0);
-  PowerMethod<CoarseVector> cPM;
-  RealD lambda_max = cPM(MrhsCoarseOp, pm_src);
-  // Chebyshev filter window [lo, hi]:
-  //   lo must sit in the spectral gap between the Nstop-th and (Nstop+1)-th
-  //   coarse eigenvalues so that only the target modes receive cosh amplification.
-  //
-  // From a pilot run (16^4 fine, 4^4 coarse, mass=0.05, hot config):
-  //   Group 1 (near-null, 24 modes): lambda in [0.002647, 0.002746]  ~= mass^2
-  //   Spectral gap: factor 60 (lambda_24/lambda_23 = 0.165/0.00275)
-  //   Group 2 (second group): lambda in [0.165, 0.179]
-  //
-  // lo = 0.02 sits in the spectral gap (factor 7x above lambda_23=0.00275,
-  // factor 8x below lambda_24=0.165).
-  //   hi = lambda_max_coarse * 1.1 ~= 2.121
-  //   y(lambda_0=0.002647)  ~ -1.016 -> T_70 ~ 1.7e5  (cosh(70*0.182))
-  //   y(lambda_23=0.002746) ~ -1.015 -> T_70 ~ 1.6e5
-  //   Relative spread across near-null cluster: ~4.3%
-  //   y(lambda_24=0.165)    ~ -0.862 -> inside [lo,hi] -> |T_70| <= 1
-  //
-  // order=71 (degree 70) is needed to give ~4% relative spread across the
-  // near-null cluster of 24 nearly-degenerate eigenvalues; order=31 (tried)
-  // gave only ~1.7% spread, insufficient for Nk=24/Nm=48 to converge.
-  // Absolute amplification ~1e5; what matters for IRL convergence is the
-  // relative spread, not the absolute value.
-  // lo=0.005 failed (T_70~53, 0/24 modes in 10 restarts).
-  // lo=0.01 worked but needed 2 restarts (13/24 then 24/24); lo=0.02 converges in 1.
-  RealD lambda_lo  = 0.02;
-
-  std::cout << GridLogMessage << "Chebyshev filter: lo=" << lambda_lo
-            << " hi=" << lambda_max*1.1 << " order=71" << std::endl;
-
-  Chebyshev<CoarseVector> IRLCheby(lambda_lo, lambda_max * 1.1, 71);
-
-  // 24 near-null modes (eigenvalues ~mass^2) converge to resid^2~1e-28
-  // in the first Lanczos restart.  The remaining modes (~0.165) are a
-  // second spectral group that needs more Krylov vectors; handle them
-  // separately once the basic HDCG solve is validated.
-  int Nk    = 24;
-  int Nm    = 48;
-  int Nstop = Nk;
-
-  GridParallelRNG CRNG(Coarse4d); CRNG.SeedFixedIntegers({13,14,15,16});
-
-  ImplicitlyRestartedBlockLanczosCoarse<CoarseVector>
-    IRL(MrhsCoarseOp, Coarse4d, CoarseMrhs, nrhs, IRLCheby,
-        Nstop, /*conv_test_interval*/1, nrhs, Nk, Nm, 1.0e-5, 10);
-
-  int Nconv;
-  std::vector<RealD>        eval(Nm);
-  std::vector<CoarseVector> evec(Nm,   Coarse4d);  // evec on f_grid (single-RHS coarse)
-  std::vector<CoarseVector> c_srcs(nrhs, Coarse4d); // src on same grid as evec
-  for (int r = 0; r < nrhs; r++) random(CRNG, c_srcs[r]);
-  IRL.calc(eval, evec, c_srcs, Nconv, LanczosType::irbl);
-
-  //--------------------------------------------------------------------
-  // HDCG solver assembly
-  //--------------------------------------------------------------------
-  MultiRHSDeflation<CoarseVector> MrhsGuesser;
-  MrhsGuesser.ImportEigenBasis(evec, eval);
-
-  // MrhsProjector maps between fine (UrbGrid) and coarse (Coarse4d) spaces
-  MultiRHSBlockProject<LatticeStaggeredFermionD> MrhsProjector;
-  MrhsProjector.Allocate(nbasis, UrbGrid, Coarse4d);
-  MrhsProjector.ImportBasis(Aggregates.subspace);
-
-  ConjugateGradient<CoarseVector> CoarseCG(5.0e-2, 5000, false);
-  DoNothingGuesser<CoarseVector>  DoNothing;
-  HPDSolver<CoarseVector> HPDSolve(MrhsCoarseOp, CoarseCG, DoNothing);
-
-  // Spectral radius of the fine operator, needed for the smoother shift.
-  // Use a random checkerboard vector (UrbGrid) as starting guess for PowerMethod.
-  LatticeStaggeredFermionD fine_pm_src(UrbGrid);
-  random(RNGrb, fine_pm_src);
-  PowerMethod<LatticeStaggeredFermionD> finePM;
-  RealD fine_lambda_max = finePM(HermOp, fine_pm_src);
-
-  // Shifted smoother: CG on (HermOp + shift*I) with shift = lambda_max / 100.
-  //
-  // The O(8) CG polynomial has 8 roots. With this shift all 8 roots lie in the
-  // interval [shift, lambda_max + shift] ~ [0.046, 4.65], so the polynomial
-  // focuses entirely on the HIGH-frequency part of the spectrum and leaves
-  // near-null modes (lambda << shift) essentially untouched (polynomial ~ 1 there).
-  //
-  // This is the right target because the coarse-grid correction always introduces
-  // high-frequency spectral leakage: the blocked coarse-grid degrees of freedom
-  // are piecewise constant across coarse cells and therefore have sharp edges at
-  // cell boundaries (like lego-block edges). Smoothness is measured by the
-  // covariant Dirac derivative, so promoting the coarse solution back to the
-  // fine grid inevitably excites high-frequency components — just as a step
-  // function always carries high-frequency Fourier content.  The smoother must
-  // repair exactly these high modes.
-  //
-  // The smoother and the coarse-grid correction are applied alternately: together
-  // they both lift the low eigenvalues and pull down the upper eigenvalues of the
-  // composite preconditioned operator, reducing the condition number seen by the
-  // outer HDCG iterations.
-  //
-  // DWF HDCG convention; using mass^2 = 0.0025 was far too small: it scattered
-  // the 8 roots over [0.005, 4.6] and diluted their effect on the high modes.
-  RealD smootherShift = fine_lambda_max / 200.0;
-  std::cout << GridLogMessage << "Smoother shift: lambda_max_fine/200 = "
-            << fine_lambda_max << "/200 = " << smootherShift << std::endl;
-  ShiftedHermOpLinearOperator<LatticeStaggeredFermionD> ShiftedOp(HermOp, smootherShift);
-  CGSmoother<LatticeStaggeredFermionD> smoother(8, ShiftedOp);
-
-  TwoLevelADEF2mrhs<LatticeStaggeredFermionD, CoarseVector>
-    HDCG(1.0e-8, 500,
-         HermOp,
-         smoother,
-         HPDSolve,   // M1 (coarse correction)
-         HPDSolve,   // Vstart (initial guess projection)
-         MrhsProjector,
-         MrhsGuesser,
-         CoarseMrhs);
-
-  //--------------------------------------------------------------------
-  // Solve: nrhs right-hand sides simultaneously
-  //--------------------------------------------------------------------
-  std::vector<LatticeStaggeredFermionD> src(nrhs, UrbGrid);
-  std::vector<LatticeStaggeredFermionD> sol(nrhs, UrbGrid);
-
-  GridParallelRNG RNGrb2(UGrid); RNGrb2.SeedFixedIntegers({17,18,19,20}); // must use full grid, not UrbGrid
-  for (int r = 0; r < nrhs; r++) {
-    random(RNGrb2, src[r]);
-    sol[r] = Zero();
-  }
-
-  //--------------------------------------------------------------------
-  // Baseline: standard single-RHS CG on HermOp (no preconditioning)
-  // Run before HDCG to establish the unpreconditioned iteration count
-  // and wall-clock time for direct comparison.
-  //--------------------------------------------------------------------
-  {
-    ConjugateGradient<LatticeStaggeredFermionD> CG(1.0e-8, 50000, false);
-    std::vector<LatticeStaggeredFermionD> cg_sol(nrhs, UrbGrid);
-    for (int r = 0; r < nrhs; r++) cg_sol[r] = Zero();
-
-    RealD t0 = usecond();
-    int total_iters = 0;
-    for (int r = 0; r < nrhs; r++) {
-      std::cout << GridLogMessage << "====== CG baseline RHS " << r
-                << " ======" << std::endl;
-      CG(HermOp, src[r], cg_sol[r]);
-      total_iters += CG.IterationsToComplete;
-    }
-    RealD t1 = usecond();
-    std::cout << GridLogMessage << "CG baseline: " << nrhs << " RHS, "
-              << total_iters << " total iterations, "
-              << (t1 - t0) / 1.0e6 << " s total, "
-              << (t1 - t0) / 1.0e6 / nrhs << " s/RHS" << std::endl;
-  }
-
-  //--------------------------------------------------------------------
-  // HDCG solve
-  //--------------------------------------------------------------------
-  HDCG(src, sol);
-
-  Grid_finalize();
-  return 0;
-}
@@ -1,373 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid
-
-    Source file: tests/debug/Test_staggered_hdcg.cc
-
-    Authors: Thomas Blum, Peter Boyle
-
-    HDCG (Hierarchical Deflation Conjugate Gradient) multigrid solver
-    for naive staggered fermions, based on arXiv:2409.03904.
-
-    Adapts the DWF HDCG infrastructure (Test_general_coarse_hdcg_phys48.cc) to:
-      - NaiveStaggeredFermion (nearest-neighbour only, no Naik 3-hop term)
-      - 4D SchurStaggeredOperator:  Mpc = m^2 - D_oe * D_eo  (hermitian, positive-definite)
-      - vColourVector fine field type (staggered has colour but no spin)
-      - NextToNearestStencilGeometry4D: 33-point coarse stencil
-
-    Stencil count: D_oe*D_eo has 2-hop fine range.  With blocking B >= 2 the coarse
-    shifts have L1-distance <= 2, giving 33 stencil points in 4D:
-      1 (identity) + 8 (+-e_mu) + 24 (+-e_mu +- e_nu).
-    NaiveStaggeredFermion has no Naik term, so any B >= 2 suffices.
-    To extend to ImprovedStaggeredFermion later, use B >= 6.
-
-    Reference: arXiv:2409.03904 (mrhs hermitian multigrid for DWF).
-
-    Usage (after build):
-      ./Test_staggered_hdcg --grid 16.16.16.16 --mpi 1.1.1.1
-
-*************************************************************************************/
-#include <Grid/Grid.h>
-#include <Grid/algorithms/iterative/AdefMrhs.h>
-#include <Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h>
-#include <Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczosCoarse.h>
-
-using namespace Grid;
-
-// Non-converging CG used as a smoother (fixed number of iterations)
-template<class Field>
-class CGSmoother : public LinearFunction<Field>
-{
-public:
-  typedef LinearOperatorBase<Field> FineOperator;
-  FineOperator &_Op;
-  int iters;
-  CGSmoother(int _iters, FineOperator &Op) : _Op(Op), iters(_iters) {}
-  void operator()(const Field &in, Field &out)
-  {
-    ConjugateGradient<Field> CG(0.0, iters, false);
-    out = Zero();
-    CG(_Op, in, out);
-  }
-};
-
-int main(int argc, char **argv)
-{
-  fprintf(stderr, "TRACE: entering main\n"); fflush(stderr);
-  Grid_init(&argc, &argv);
-  fprintf(stderr, "TRACE: Grid_init done\n"); fflush(stderr);
-
-  //--------------------------------------------------------------------
-  // Parameters — tune for production
-  //--------------------------------------------------------------------
-  const int nbasis = 24;   // near-null space dimension
-  const int cb     = 0;    // even checkerboard
-
-  RealD mass = 0.00184;
-
-  // NaiveStaggeredFermion: nearest-neighbour hop only (no Naik term).
-  // c1 = coefficient of the hopping term (1.0 = standard normalisation).
-  // u0 = tadpole factor (1.0 = no tadpole improvement).
-  RealD c1 = 1.0;
-  RealD u0 = 1.0;
-
-  //--------------------------------------------------------------------
-  // Grids
-  // Fine:   UGrid (4D full), UrbGrid (4D red-black)
-  // Coarse: Coarse4d  with dimensions = GridDefaultLatt() / Block
-  //
-  // Recommended: GridDefaultLatt() >= 16^4, Block = {4,4,4,4}
-  // NaiveStaggeredFermion works with any Block >= {2,2,2,2}
-  //--------------------------------------------------------------------
-  fprintf(stderr, "TRACE: making UGrid\n"); fflush(stderr);
-  GridCartesian *UGrid = SpaceTimeGrid::makeFourDimGrid(
-      GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi());
-  fprintf(stderr, "TRACE: making UrbGrid\n"); fflush(stderr);
-  GridRedBlackCartesian *UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
-
-  Coordinate Block({4, 4, 4, 4});
-  Coordinate clatt = GridDefaultLatt();
-  for (int d = 0; d < clatt.size(); d++) clatt[d] /= Block[d];
-  Coordinate csimd = GridDefaultSimd(Nd, vComplex::Nsimd());
-  Coordinate cmpi  = GridDefaultMpi();
-  fprintf(stderr, "TRACE: making Coarse4d clatt=%d %d %d %d simd=%d %d %d %d mpi=%d %d %d %d Nsimd=%d\n",
-          clatt[0],clatt[1],clatt[2],clatt[3],
-          csimd[0],csimd[1],csimd[2],csimd[3],
-          cmpi[0],cmpi[1],cmpi[2],cmpi[3],
-          (int)vComplex::Nsimd()); fflush(stderr);
-
-  GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, csimd, cmpi);
-  fprintf(stderr, "TRACE: Coarse4d made\n"); fflush(stderr);
-
-  //--------------------------------------------------------------------
-  // RNG + gauge field
-  //--------------------------------------------------------------------
-  fprintf(stderr, "TRACE: RNG4\n"); fflush(stderr);
-  GridParallelRNG RNG4(UGrid);    RNG4.SeedFixedIntegers({1,2,3,4});
-  fprintf(stderr, "TRACE: RNGrb\n"); fflush(stderr);
-  GridParallelRNG RNGrb(UGrid); RNGrb.SeedFixedIntegers({5,6,7,8}); // must use full grid, not UrbGrid
-  fprintf(stderr, "TRACE: Umu\n"); fflush(stderr);
-  LatticeGaugeField Umu(UGrid);
-  int HotStart = 0;
-  if ( HotStart ) {
-    fprintf(stderr, "TRACE: HotConfig\n"); fflush(stderr);
-    SU<Nc>::HotConfiguration(RNG4, Umu);
-  } else { 
-    FieldMetaData header;
-    std::string file("./configuration.ildg");
-    IldgReader IR;
-    IR.open(file);
-    IR.readConfiguration(Umu,header);
-    IR.close();
-  }
-  
-  fprintf(stderr, "TRACE: NaiveStaggeredFermionD\n"); fflush(stderr);
-  NaiveStaggeredFermionD Ds(Umu, *UGrid, *UrbGrid, mass, c1, u0);
-  fprintf(stderr, "TRACE: SchurStaggeredOperator\n"); fflush(stderr);
-  SchurStaggeredOperator<NaiveStaggeredFermionD, LatticeStaggeredFermionD> HermOp(Ds);
-  fprintf(stderr, "TRACE: HermOp done\n"); fflush(stderr);
-
-  //--------------------------------------------------------------------
-  // Subspace: inverse-iteration near-null vectors
-  //
-  // CreateSubspace applies CG (4 solves, tol=1e-4) to random noise vectors,
-  // converging naturally to the low modes of HermOp without needing spectral
-  // bound tuning.  Switch to CreateSubspaceChebyshevNew once the spectrum is
-  // well characterised (hi ~ 5.0 for naive staggered SchurStaggeredOperator).
-  //--------------------------------------------------------------------
-  typedef Aggregation<vColourVector, vTComplex, nbasis> Subspace;
-  Subspace Aggregates(Coarse4d, UrbGrid, cb);
-
-  Aggregates.CreateSubspace(RNGrb, HermOp);
-  Aggregates.Orthogonalise();
-
-  //--------------------------------------------------------------------
-  // Coarse geometry: NextToNearestStencilGeometry4D
-  //   hops=2  ->  33 stencil points in 4D
-  //--------------------------------------------------------------------
-  NextToNearestStencilGeometry4D geom(Coarse4d);
-
-  std::cout << GridLogMessage << "Coarse stencil: " << geom.npoint << " points" << std::endl;
-
-  //--------------------------------------------------------------------
-  // Single-RHS coarse operator (used for correctness check below)
-  //--------------------------------------------------------------------
-  typedef GeneralCoarsenedMatrix<vColourVector, vTComplex, nbasis> LittleDiracOp;
-  typedef LittleDiracOp::CoarseVector CoarseVector;
-
-  LittleDiracOp LDO(geom, UrbGrid, Coarse4d);
-  LDO.CoarsenOperator(HermOp, Aggregates);
-
-  //--------------------------------------------------------------------
-  // Correctness check: P M_fine P^T c  ≈  M_coarse c
-  //
-  // Promote a random coarse vector into the fine subspace, apply the
-  // fine operator, project back, and compare with the coarse operator
-  // applied directly.  Error should be at the level of subspace
-  // approximation quality (smaller = better basis vectors).
-  //--------------------------------------------------------------------
-  {
-    GridParallelRNG RNGc(Coarse4d); RNGc.SeedFixedIntegers({9,10,11,12});
-    CoarseVector c_src(Coarse4d), c_ldop(Coarse4d), c_proj(Coarse4d);
-    random(RNGc, c_src);
-
-    LatticeStaggeredFermionD f_v(UrbGrid), f_Mv(UrbGrid);
-    Aggregates.PromoteFromSubspace(c_src, f_v);
-    HermOp.Op(f_v, f_Mv);
-    Aggregates.ProjectToSubspace(c_proj, f_Mv);
-
-    LDO.M(c_src, c_ldop);
-
-    c_proj -= c_ldop;
-    RealD err = norm2(c_proj) / norm2(c_ldop);
-    std::cout << GridLogMessage
-              << "Coarsen check |P*M_fine - M_coarse| / |M_coarse| = " << err << std::endl;
-  }
-
-  //--------------------------------------------------------------------
-  // Multi-RHS coarse grid
-  //
-  // The extra leading dimension holds nrhs right-hand sides packed into
-  // SIMD lanes, matching the pattern of Test_general_coarse_hdcg_phys48.
-  //--------------------------------------------------------------------
-  const int nrhs = vComplex::Nsimd() * 2;
-
-  Coordinate mpi   = GridDefaultMpi();
-  Coordinate rhMpi ({1,    mpi[0],  mpi[1],  mpi[2],  mpi[3]});
-  Coordinate rhLatt({nrhs, clatt[0], clatt[1], clatt[2], clatt[3]});
-  Coordinate rhSimd({vComplex::Nsimd(), 1, 1, 1, 1});
-
-  GridCartesian *CoarseMrhs = new GridCartesian(rhLatt, rhSimd, rhMpi);
-
-  typedef MultiGeneralCoarsenedMatrix<vColourVector, vTComplex, nbasis> MultiCoarseOp;
-  MultiCoarseOp mrhs(geom, CoarseMrhs);
-  mrhs.CoarsenOperator(HermOp, Aggregates, Coarse4d);
-
-  //--------------------------------------------------------------------
-  // Coarse-grid Lanczos for deflation
-  //--------------------------------------------------------------------
-  typedef HermitianLinearOperator<MultiCoarseOp, CoarseVector> MrhsHermOp;
-  MrhsHermOp MrhsCoarseOp(mrhs);
-
-  // Estimate spectral bounds for Lanczos Chebyshev filter
-  CoarseVector pm_src(CoarseMrhs); pm_src = ComplexD(1.0);
-  PowerMethod<CoarseVector> cPM;
-  RealD lambda_max = cPM(MrhsCoarseOp, pm_src);
-  // Chebyshev filter window [lo, hi]:
-  //   lo must sit in the spectral gap between the Nstop-th and (Nstop+1)-th
-  //   coarse eigenvalues so that only the target modes receive cosh amplification.
-  //
-  // From a pilot run (16^4 fine, 4^4 coarse, mass=0.05, hot config):
-  //   Group 1 (near-null, 24 modes): lambda in [0.002647, 0.002746]  ~= mass^2
-  //   Spectral gap: factor 60 (lambda_24/lambda_23 = 0.165/0.00275)
-  //   Group 2 (second group): lambda in [0.165, 0.179]
-  //
-  // lo = 0.02 sits in the spectral gap (factor 7x above lambda_23=0.00275,
-  // factor 8x below lambda_24=0.165).
-  //   hi = lambda_max_coarse * 1.1 ~= 2.121
-  //   y(lambda_0=0.002647)  ~ -1.016 -> T_70 ~ 1.7e5  (cosh(70*0.182))
-  //   y(lambda_23=0.002746) ~ -1.015 -> T_70 ~ 1.6e5
-  //   Relative spread across near-null cluster: ~4.3%
-  //   y(lambda_24=0.165)    ~ -0.862 -> inside [lo,hi] -> |T_70| <= 1
-  //
-  // order=71 (degree 70) is needed to give ~4% relative spread across the
-  // near-null cluster of 24 nearly-degenerate eigenvalues; order=31 (tried)
-  // gave only ~1.7% spread, insufficient for Nk=24/Nm=48 to converge.
-  // Absolute amplification ~1e5; what matters for IRL convergence is the
-  // relative spread, not the absolute value.
-  // lo=0.005 failed (T_70~53, 0/24 modes in 10 restarts).
-  // lo=0.01 worked but needed 2 restarts (13/24 then 24/24); lo=0.02 converges in 1.
-  RealD lambda_lo  = 0.02;
-
-  std::cout << GridLogMessage << "Chebyshev filter: lo=" << lambda_lo
-            << " hi=" << lambda_max*1.1 << " order=71" << std::endl;
-
-  Chebyshev<CoarseVector> IRLCheby(lambda_lo, lambda_max * 1.1, 71);
-
-  // 24 near-null modes (eigenvalues ~mass^2) converge to resid^2~1e-28
-  // in the first Lanczos restart.  The remaining modes (~0.165) are a
-  // second spectral group that needs more Krylov vectors; handle them
-  // separately once the basic HDCG solve is validated.
-  int Nk    = 24;
-  int Nm    = 48;
-  int Nstop = Nk;
-
-  GridParallelRNG CRNG(Coarse4d); CRNG.SeedFixedIntegers({13,14,15,16});
-
-  ImplicitlyRestartedBlockLanczosCoarse<CoarseVector>
-    IRL(MrhsCoarseOp, Coarse4d, CoarseMrhs, nrhs, IRLCheby,
-        Nstop, /*conv_test_interval*/1, nrhs, Nk, Nm, 1.0e-5, 10);
-
-  int Nconv;
-  std::vector<RealD>        eval(Nm);
-  std::vector<CoarseVector> evec(Nm,   Coarse4d);  // evec on f_grid (single-RHS coarse)
-  std::vector<CoarseVector> c_srcs(nrhs, Coarse4d); // src on same grid as evec
-  for (int r = 0; r < nrhs; r++) random(CRNG, c_srcs[r]);
-  IRL.calc(eval, evec, c_srcs, Nconv, LanczosType::irbl);
-
-  //--------------------------------------------------------------------
-  // HDCG solver assembly
-  //--------------------------------------------------------------------
-  MultiRHSDeflation<CoarseVector> MrhsGuesser;
-  MrhsGuesser.ImportEigenBasis(evec, eval);
-
-  // MrhsProjector maps between fine (UrbGrid) and coarse (Coarse4d) spaces
-  MultiRHSBlockProject<LatticeStaggeredFermionD> MrhsProjector;
-  MrhsProjector.Allocate(nbasis, UrbGrid, Coarse4d);
-  MrhsProjector.ImportBasis(Aggregates.subspace);
-
-  ConjugateGradient<CoarseVector> CoarseCG(5.0e-2, 5000, false);
-  DoNothingGuesser<CoarseVector>  DoNothing;
-  HPDSolver<CoarseVector> HPDSolve(MrhsCoarseOp, CoarseCG, DoNothing);
-
-  // Spectral radius of the fine operator, needed for the smoother shift.
-  // Use a random checkerboard vector (UrbGrid) as starting guess for PowerMethod.
-  LatticeStaggeredFermionD fine_pm_src(UrbGrid);
-  random(RNGrb, fine_pm_src);
-  PowerMethod<LatticeStaggeredFermionD> finePM;
-  RealD fine_lambda_max = finePM(HermOp, fine_pm_src);
-
-  // Shifted smoother: CG on (HermOp + shift*I) with shift = lambda_max / 100.
-  //
-  // The O(8) CG polynomial has 8 roots. With this shift all 8 roots lie in the
-  // interval [shift, lambda_max + shift] ~ [0.046, 4.65], so the polynomial
-  // focuses entirely on the HIGH-frequency part of the spectrum and leaves
-  // near-null modes (lambda << shift) essentially untouched (polynomial ~ 1 there).
-  //
-  // This is the right target because the coarse-grid correction always introduces
-  // high-frequency spectral leakage: the blocked coarse-grid degrees of freedom
-  // are piecewise constant across coarse cells and therefore have sharp edges at
-  // cell boundaries (like lego-block edges). Smoothness is measured by the
-  // covariant Dirac derivative, so promoting the coarse solution back to the
-  // fine grid inevitably excites high-frequency components — just as a step
-  // function always carries high-frequency Fourier content.  The smoother must
-  // repair exactly these high modes.
-  //
-  // The smoother and the coarse-grid correction are applied alternately: together
-  // they both lift the low eigenvalues and pull down the upper eigenvalues of the
-  // composite preconditioned operator, reducing the condition number seen by the
-  // outer HDCG iterations.
-  //
-  // DWF HDCG convention; using mass^2 = 0.0025 was far too small: it scattered
-  // the 8 roots over [0.005, 4.6] and diluted their effect on the high modes.
-  RealD smootherShift = fine_lambda_max / 200.0;
-  std::cout << GridLogMessage << "Smoother shift: lambda_max_fine/200 = "
-            << fine_lambda_max << "/200 = " << smootherShift << std::endl;
-  ShiftedHermOpLinearOperator<LatticeStaggeredFermionD> ShiftedOp(HermOp, smootherShift);
-  CGSmoother<LatticeStaggeredFermionD> smoother(8, ShiftedOp);
-
-  TwoLevelADEF2mrhs<LatticeStaggeredFermionD, CoarseVector>
-    HDCG(1.0e-8, 500,
-         HermOp,
-         smoother,
-         HPDSolve,   // M1 (coarse correction)
-         HPDSolve,   // Vstart (initial guess projection)
-         MrhsProjector,
-         MrhsGuesser,
-         CoarseMrhs);
-
-  //--------------------------------------------------------------------
-  // Solve: nrhs right-hand sides simultaneously
-  //--------------------------------------------------------------------
-  std::vector<LatticeStaggeredFermionD> src(nrhs, UrbGrid);
-  std::vector<LatticeStaggeredFermionD> sol(nrhs, UrbGrid);
-
-  GridParallelRNG RNGrb2(UGrid); RNGrb2.SeedFixedIntegers({17,18,19,20}); // must use full grid, not UrbGrid
-  for (int r = 0; r < nrhs; r++) {
-    random(RNGrb2, src[r]);
-    sol[r] = Zero();
-  }
-
-  //--------------------------------------------------------------------
-  // Baseline: standard single-RHS CG on HermOp (no preconditioning)
-  // Run before HDCG to establish the unpreconditioned iteration count
-  // and wall-clock time for direct comparison.
-  //--------------------------------------------------------------------
-  {
-    ConjugateGradient<LatticeStaggeredFermionD> CG(1.0e-8, 100000, false);
-    std::vector<LatticeStaggeredFermionD> cg_sol(nrhs, UrbGrid);
-    for (int r = 0; r < nrhs; r++) cg_sol[r] = Zero();
-
-    RealD t0 = usecond();
-    int total_iters = 0;
-    for (int r = 0; r < nrhs; r++) {
-      std::cout << GridLogMessage << "====== CG baseline RHS " << r
-                << " ======" << std::endl;
-      CG(HermOp, src[r], cg_sol[r]);
-      total_iters += CG.IterationsToComplete;
-    }
-    RealD t1 = usecond();
-    std::cout << GridLogMessage << "CG baseline: " << nrhs << " RHS, "
-              << total_iters << " total iterations, "
-              << (t1 - t0) / 1.0e6 << " s total, "
-              << (t1 - t0) / 1.0e6 / nrhs << " s/RHS" << std::endl;
-  }
-
-  //--------------------------------------------------------------------
-  // HDCG solve
-  //--------------------------------------------------------------------
-  HDCG(src, sol);
-
-  Grid_finalize();
-  return 0;
-}
Author	SHA1	Message	Date
Peter Boyle	1f55c13bc4	Test script for Frontier on I/O	2026-08-11 13:37:50 -04:00
Peter Boyle	d68b111d06	Switch on the faster IO be default	2026-08-11 13:36:50 -04:00
Peter Boyle	41f5a02204	IO updates for AllToAllV aggregation of contiguous chunks. May finally clean up the poor MPI2 IO performance issue that has been persistent.	2026-08-11 13:33:04 -04:00
Peter Boyle	7647576863	Sourceme's for frontier	2026-06-26 11:50:44 -04:00
Peter Boyle	f11ba18df2	Remove accelerator_inline on CPU only code	2026-06-17 20:47:15 +02:00
Peter Boyle	cf8587e401	Update booster compiule	2026-06-17 20:46:53 +02:00
Peter Boyle	7dd35ef749	Make disable accelerator aware mpi compile for CUDA	2026-06-17 20:46:14 +02:00
Peter Boyle	41e570ddce	Annoying old CPU perfmon code should be removed or deprecated as not worth maintaining	2026-06-17 20:45:32 +02:00
Peter Boyle	a452131b50	Print improvement	2026-06-17 20:45:04 +02:00
Peter Boyle	4e49ca55ab	Merge branch 'develop' of https://github.com/paboyle/Grid into develop	2026-06-16 11:20:12 -07:00
Peter Boyle	c3f4474401	Adding Mattia's memory leak test	2026-06-16 11:19:36 -07:00
Peter Boyle	3d3eff86f3	Modify move assignment operator to be noexcept Add noexcept specifier to move assignment operator.	2026-06-11 09:44:24 -04:00
Peter Boyle	fc9f154ac1	Modify Lattice move constructor to be noexcept Add noexcept specifier to move constructor for Lattice class.	2026-06-11 09:40:06 -04:00
Peter Boyle	4aa0bca4dc	Change sum operation to use gpucub mistake in PR from Chris Updated the sum operation definition for GPU reduction to use gpucub instead of cub.	2026-06-01 14:12:25 -04:00
Peter Boyle	905da6f083	Merge branch 'feature/reduction-reorganisation' into develop	2026-05-27 21:01:30 -04:00
Peter Boyle	12e3499b6d	Updated rocm 7 compile for ORNL	2026-05-21 12:28:42 -04:00
Peter Boyle	9576011011	Changed setup for ROCM 7, nasty LD_LIBRARY_PATH issues were committing evils	2026-05-21 12:28:04 -04:00
Peter Boyle	155b34c1aa	File list lost	2026-05-21 12:06:01 -04:00
Peter BoyleandClaude Sonnet 4.6	982ffe9ebe	Lattice_reduction_gpu: demote timing logs to Debug, disable by default skills/mpi-heterogeneous: add Bug Class 4 for Frontier GTL/libamdhip64 ABI mismatch Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-21 12:05:36 -04:00
Peter BoyleandClaude Sonnet 4.6	0251ecaeab	Test_planned_fft: fix PlannedFFT template parameter to use ::vector_object Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-20 18:13:38 -04:00
Peter BoyleandClaude Sonnet 4.6	372a27d645	tests: add Test_planned_fft exercising PlannedFFT<vobj> Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-20 17:59:24 -04:00
Peter BoyleandClaude Sonnet 4.6	72b4a061f3	tests/fft: remove PlanDestroy calls (FFT handles plans per-call) Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-20 17:54:41 -04:00
Peter BoyleandClaude Sonnet 4.6	29198efabe	FFT: add FFTbase, PlannedFFT; factor FFT_dim_execute free function Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-20 17:53:17 -04:00
Peter BoyleandClaude Sonnet 4.6	50aa51f93a	debug: add Test_hipfft_repro — reproducer for hipFFT PARSE_ERROR on ROCm 7 Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-19 22:27:27 -04:00
Peter BoyleandClaude Sonnet 4.6	79ccc81a86	tests/debug: add G=4 to hipfft fail reproducer Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-19 22:21:52 -04:00
Peter BoyleandClaude Sonnet 4.6	3f0fdbb597	tests/debug: test hipMemset variant before cache is populated Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-19 22:10:16 -04:00
Peter BoyleandClaude Sonnet 4.6	ea57bd8f03	tests/debug: extend hipfft fail reproducer with hipMemset and sync variants Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-19 22:02:02 -04:00
Peter BoyleandClaude Sonnet 4.6	bdba5b8403	FFT: use host stack buffer in PlanCreate, not deviceVector Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-19 21:49:06 -04:00
Peter BoyleandClaude Sonnet 4.6	58cc6ca9c0	tests/debug: add minimal hipfft ordering bug fail/pass pair Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-19 21:48:23 -04:00
Peter BoyleandClaude Sonnet 4.6	e5996b440d	tests/debug: test plan-before-malloc vs malloc-before-plan ordering Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-19 21:40:17 -04:00
Peter BoyleandClaude Sonnet 4.6	ad9d03fd85	tests/debug: extend hipfft reproducer with Grid-realistic howmany and exec tests Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-19 19:19:59 -04:00
Peter BoyleandClaude Sonnet 4.6	4de160ce20	tests/debug: add minimal hipfft plan-creation reproducer Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-19 17:52:59 -04:00
Peter Boyle	fc8c8ce6e7	FFT HIP: use hipfftCreate+hipfftMakePlanMany instead of hipfftPlanMany	2026-05-19 17:29:28 -04:00
Peter Boyle	ddbb7f07c8	FFT: pass nullptr for inembed/onembed in hipfftPlanMany to avoid HIPFFT_PARSE_ERROR	2026-05-19 17:15:21 -04:00
Peter BoyleandClaude Sonnet 4.6	1e29c59bcc	FFT: cache plans per vobj type across calls Plans are created lazily on the first FFT_dim call and reused for all subsequent calls on the same FFT object. PlanCreate<vobj>() can be called explicitly to pre-warm the cache. PlanDestroy() must be called before switching to a different vobj type; the destructor cleans up any live plans automatically. Update Test_fft.cc and Test_fftf.cc to call PlanDestroy() between the LatticeComplex and LatticeSpinMatrix sections that reuse the same FFT object. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-19 15:12:10 -04:00
Peter BoyleandClaude Sonnet 4.6	b6abdc3845	Accelerator: lower default accelerator_threads from 16 to 8 Benchmark_dwf_fp32 on MI250X GCD: 1.7 TF/s at nt=8, ~300 GF/s at nt=16. With Nsimd=8 (fp32, GEN_SIMD_WIDTH=64B), nt=8 gives exactly 64 threads = one full AMD wavefront. Higher values double register demand per block and hit a register-pressure cliff for stencil kernels. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-19 13:41:03 -04:00
Peter Boyle	2fadd8bb62	Accelerator: raise default accelerator_threads from 2 to 16	2026-05-19 10:15:53 -04:00
Peter BoyleandClaude Sonnet 4.6	60df2dd5d0	skills: add gpu-memory-performance.md Documents the acceleratorThreads() default=2 trap, LambdaApply thread mapping, coalescedRead/Write idiom, when to use __global__ vs accelerator_for, and fused vs staged HBM access patterns. Includes observed MI250X numbers from LatticePropagatorD reduction (50 → 297 → 546 GB/s progression). Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-19 10:03:32 -04:00
Peter BoyleandClaude Sonnet 4.6	66b529b345	sumD_gpu_reduce_words: fuse pack+reduce into single packReduceKernel Replace the two-kernel pack+reduce sequence with a single fused kernel packReduceKernel<R> that reads R words of each vobj at offset 'base' and accumulates directly into iVector<iScalar<scalarD>,R>, eliminating the intermediate bundle buffer entirely. HBM access per word-group drops from 3x (pack-read + pack-write + reduce-read) to 1x. Thread count comes from getNumBlocksAndThreads (warpSize..256) rather than acceleratorThreads(), so occupancy is correct regardless of the --accelerator-threads setting. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-19 09:46:43 -04:00
Peter Boyle	1304172a93	Modified repack	2026-05-19 08:53:13 -04:00
Peter BoyleandClaude Sonnet 4.6	1315d4604d	Enable GRID_REDUCTION_TIMING unconditionally Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-18 22:14:00 -04:00
Peter BoyleandClaude Sonnet 4.6	a31af31328	Lattice_reduction_gpu: add GRID_REDUCTION_TIMING instrumentation Uncomment #define GRID_REDUCTION_TIMING to enable per-phase timing output: sumD_gpu_reduce_words: pack time (accelerator_for) per R and base sumD_gpu_small: reduceKernel+barrier time and D2H time separately sumD_gpu_large: total wall time across all word groups This lets us identify whether the large-type bottleneck is in the pack kernel, the shared-memory reduction kernel, the barrier, or the D2H. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-18 22:13:30 -04:00
Peter BoyleandClaude Sonnet 4.6	26c3c7d8f9	sumD_gpu_large: radix-12 word-bundle reduction replacing radix-1 Replace the word-by-word loop (one kernel launch per scalar word) with sumD_gpu_reduce_words<R> which packs R consecutive vector_type words per site into iVector<iScalar<vector>,R>, then calls the existing sumD_gpu_small shared-memory kernel once for the whole bundle. Dispatch: radix-12 first, radix-4 for the remainder < 12, radix-1 for any final < 4 words. For LatticePropagator (144 words = 12x12), this reduces the kernel-launch count from 144 to 12 -- a 12x reduction. Bundle::Nsimd() inherits from vector_type so sumD_gpu_small handles SIMD lane extraction and double-precision promotion identically to the scalar word case. sizeof(Bundle::scalar_objectD) = R*16 <= 192 B; well within sharedMemPerBlock on all supported devices. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-18 21:56:45 -04:00
Peter BoyleandClaude Sonnet 4.6	0650d7c7eb	Lattice_reduction_sycl: fix double-precision accumulation in sumD_gpu_tensor Accumulate in sobjD throughout rather than accumulating in sobj and converting the final sum. For float fields this matters: summing N floats then casting loses O(N*eps_float) relative precision vs accumulating in double from the start. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-18 21:53:40 -04:00
Peter BoyleandClaude Sonnet 4.6	068f95ad2d	Revert to hand-rolled reduction; drop Lattice_reduction_gpu_cub.h Remove the CUB/hipCUB direction entirely. Restore Lattice_reduction_gpu.h, Lattice_reduction_sycl.h, and Lattice_reduction.h to the state before the CUB rewrite (commit `969b0a39`), recovering the original primary function names (sumD_gpu_small, sumD_gpu_large, sumD_gpu, sum_gpu, sum_gpu_large) and the hand-rolled shared-memory reduction kernel. Delete Lattice_reduction_gpu_cub.h. Update Test_reduction to remove the old/new comparison sections that depended on sum_gpu_old. The lesson: CUB DeviceReduce is slower than the hand-rolled kernel for small types, and the smem sizing problem for the extraction pass has no clean solution within the accelerator_for abstraction. The right improvement is a higher radix (12 then 4) in sumD_gpu_large, applied directly to the existing hand-rolled kernel. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-18 21:52:18 -04:00
Peter BoyleandClaude Sonnet 4.6	f4fbf7c9ca	sumD_gpu_direct: revert to per-lane write; CUB handles Nsimdosites inputs Benchmarking showed the shared-memory lane-summation approach (`843d6497`) was slower than writing each SIMD lane individually and letting CUB reduce the full nlanes = ositesNsimd array. CUB's device reduce is more efficient over the larger input than the smem overhead + serialised lane-0 summation. The smem approach also required overriding acceleratorThreads() to avoid the block-size sizing problem. Restore the simpler per-lane path. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-18 21:23:15 -04:00
Peter BoyleandClaude Sonnet 4.6	843d6497b2	sumD_gpu_direct: shared-memory lane reduction with acceleratorThreads(1) Set acceleratorThreads to 1 before the extraction kernel so that dim3(nsimd,1,1) blocks give exactly one site group per block and __shared__ sobjD smem[nsimd] is correctly sized without depending on the runtime acceleratorThreads() value. threadIdx.x (acceleratorSIMTlane) indexes the SIMD lane for coalesced reads; lane 0 sums smem[0..nsimd-1] and writes one sobjD per site. CUB then reduces osites elements instead of osites*nsimd, reducing both store traffic and CUB work by Nsimd. acceleratorSynchronise() (warp-level) suffices since nsimd < warpSize. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-18 21:08:10 -04:00
Peter BoyleandClaude Sonnet 4.6	747c167658	sumD_gpu_direct: one thread per SIMD lane using extractLane Replaces one thread per outer site calling Reduce() (sequential Nsimd-wide loop) with one thread per lane calling extractLane() — O(1) per thread. CUB now reduces over osites*Nsimd elements. Avoids serial lane reduction but leaves the per-lane sobjD store stride as a known remaining concern. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-18 16:21:50 -04:00
Peter BoyleandClaude Sonnet 4.6	fca2c5dba0	Lattice_reduction_gpu_cub: define GRID_REDUCTION_TIMING in header Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-18 14:54:08 -04:00
Peter BoyleandClaude Sonnet 4.6	e12bc7f07c	Lattice_reduction_gpu_cub: add GRID_REDUCTION_TIMING instrumentation Guards accelerator_for and CUB DeviceReduce calls in sumD_gpu_direct and sumD_gpu_large with #ifdef GRID_REDUCTION_TIMING to isolate where time is spent in each path. Large path accumulates across all groups and prints totals with words/nfull/rem context. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-18 14:23:44 -04:00
Peter BoyleandClaude Sonnet 4.6	dc6ae51cab	Lattice_reduction_gpu_cub: replace WordBundle4 with iVector<iScalar<scalarD>,4> WordBundle4 was redundant with Grid's existing tensor infrastructure. iVector<iScalar<scalarD>,4> already provides accelerator_inline operator+, zeroit(), and sycl::is_device_copyable — no new type needed. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-18 13:55:28 -04:00
Peter BoyleandClaude Sonnet 4.6	baa70d8ec9	Test_reduction: add timing benchmark for new vs old reduction paths Reports us/call and GB/s for sum_gpu (CUB/sycl::reduction) and sum_gpu_old (hand-rolled shared-memory) for each field type, with 5-call warmup and 100-call timed loop. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-18 12:31:13 -04:00
Peter BoyleandClaude Sonnet 4.6	c93b338bdd	skills: HPC battle-hardening skill files for GPU+MPI correctness Six skill files encoding expertise for making codebases robust on problematic HPC systems, covering: correctness verification (double-run, fingerprinting, flight recorder), hang diagnosis, GPU runtime correctness (premature barrier, infinite poll), MPI correctness on heterogeneous systems (device buffer aliasing, AARCH64 PLT corruption, deterministic reductions), compiler validation, and communication/computation overlap pipeline design. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-18 12:10:44 -04:00
Peter BoyleandClaude Sonnet 4.6	c0472aa0ec	Test_reduction: use separate float and double grids Float fields require a grid constructed with vComplexF::Nsimd(); using a double grid causes grid->_gsites to undercount the sites in float vobjF, making the constant-field expected value wrong. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-18 12:09:35 -04:00
Peter BoyleandClaude Sonnet 4.6	09552cfd73	Rename scalarNorm2 to squaredSum in Test_reduction.cc The function computes \|sum\|^2 — the squared magnitude of an aggregate sum — not a norm. squaredSum makes clear that squaring is applied to the sum, not to individual site values before summing, distinguishing it from sumOfSquares (the squared L2 norm). Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-15 23:15:11 -04:00
Peter BoyleandClaude Sonnet 4.6	003fec509c	Fix Zero() used on thrust::complex in WordBundle4 initialisation Grid's Zero() sentinel is not assignable to thrust::complex<double>; use scalarD(0) instead. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-15 18:10:17 -04:00
Peter BoyleandClaude Sonnet 4.6	773a82d87f	Reinstate large/small dispatch in CUB reduction path; radix-4 word-bundle for large types rocPRIM's DeviceReduce requires warpSize(64) threads each holding one element in shared memory, so sizeof(T)64 must fit in sharedMemPerBlock. LatticePropagator::scalar_objectD is 2304 bytes (642304 = 147 KB), exceeding the budget and triggering a compile-time static_assert in limit_block_size. Introduce sumD_gpu_direct (the original direct-CUB path, safe for small types) and a new sumD_gpu_large that groups the vobj's vector_type words in bundles of 4, reducing each bundle as WordBundle4<scalarD> (64 bytes, 64*64 = 4 KB — always within budget). If words % 4 != 0, the final partial bundle is zero-padded. sumD_gpu dispatches at compile time via if constexpr on sizeof(sobjD) > 512. For LatticePropagator (144 words) this gives 36 CUB launches instead of 144. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-15 16:55:58 -04:00
Peter BoyleandClaude Sonnet 4.6	286c29d6fb	Add Test_reduction to tests/debug Tests the new CUB/hipCUB/SYCL lattice reduction (sum_gpu) against the preserved hand-rolled implementation (sum_gpu_old) for LatticeComplexF/D, LatticeColourMatrixF/D and LatticePropagatorF/D. Part a) gaussian random field: checks that old and new agree to within float/double roundoff tolerance. Part b) constant field (= 1.0, identity-matrix init): verifies innerProduct(sum, sum) = Ncomp * V^2 where Ncomp counts the nonzero diagonal scalar components per site (1 / Nc / Ns*Nc respectively). Make.inc is auto-generated by scripts/filelist on bootstrap and is not tracked; the new .cc file is all that is needed. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-15 14:31:33 -04:00
Peter BoyleandClaude Sonnet 4.6	969b0a3922	Rewrite lattice GPU reduction to use CUB, hipCUB, and SYCL reduction Replace hand-rolled shared-memory reduction kernels (reduceBlock/reduceBlocks/ reduceKernel) and the global device variable retirementCount with a unified CUB/hipCUB DeviceReduce::Reduce path for CUDA/HIP and sycl::reduction for SYCL. No small/large split is needed: both CUB and sycl::reduction handle arbitrary object sizes internally. Old implementations preserved as sum_gpu_old / sumD_gpu_old etc. in the original files for regression testing on GPU hardware. Also add CLAUDE.md with build, test, and architecture guidance. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>	2026-05-15 13:41:56 -04:00