IO updates for AllToAllV aggregation of contiguous chunks.

May finally clean up the poor MPI2 IO performance issue that has been persistent.

CLAUDE.md

@@ -39,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)
@@ -75,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
 
@@ -85,14 +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.
+- 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`).

Grid/communicator/Communicator_base.h

@@ -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)
   {

Grid/communicator/Communicator_mpi3.cc

@@ -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);

Grid/communicator/Communicator_none.cc

@@ -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){}

Grid/parallelIO/BinaryIO.cc

@@ -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;

Grid/parallelIO/BinaryIO.h

@@ -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);

tests/IO/Test_aggregate_io.cc

@@ -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;
+}