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>

CLAUDE.md

@@ -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`
@@ -99,17 +96,3 @@ GPU support is injected via macros (`accelerator_for`, `accelerator_for2dNB`). T
 - 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/`.
 - 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 |

Grid/algorithms/FFT.h

@@ -1,6 +1,6 @@
 /*************************************************************************************
 
-    Grid physics library, www.github.com/paboyle/Grid
+    Grid physics library, www.github.com/paboyle/Grid 
 
     Source file: ./lib/Cshift.h
 
@@ -65,16 +65,17 @@ public:
   typedef hipfftDoubleComplex FFTW_scalar;
   typedef hipfftHandle        FFTW_plan;
   static FFTW_plan fftw_plan_many_dft(int rank, int *n,int howmany,
-				      FFTW_scalar *in, int *inembed,
-				      int istride, int idist,
-				      FFTW_scalar *out, int *onembed,
-				      int ostride, int odist,
+				      FFTW_scalar *in, int *inembed,		
+				      int istride, int idist,		
+				      FFTW_scalar *out, int *onembed,		
+				      int ostride, int odist,		
 				      int sign, unsigned flags) {
     FFTW_plan p;
     auto rv = hipfftPlanMany(&p,rank,n,n,istride,idist,n,ostride,odist,HIPFFT_Z2Z,howmany);
     GRID_ASSERT(rv==HIPFFT_SUCCESS);
     return p;
-  }
+  }	  
+    
   inline static void fftw_execute_dft(const FFTW_plan p,FFTW_scalar *in,FFTW_scalar *out, int sign) {
     hipfftResult rv;
     if ( sign == forward ) rv =hipfftExecZ2Z(p,in,out,HIPFFT_FORWARD);
@@ -82,25 +83,29 @@ public:
     accelerator_barrier();
     GRID_ASSERT(rv==HIPFFT_SUCCESS);
   }
-  inline static void fftw_destroy_plan(const FFTW_plan p) { hipfftDestroy(p); }
+  inline static void fftw_destroy_plan(const FFTW_plan p) {
+    hipfftDestroy(p);
+  }
 };
 template<> struct FFTW<ComplexF> {
 public:
   static const int forward=FFTW_FORWARD;
   static const int backward=FFTW_BACKWARD;
-  typedef hipfftComplex FFTW_scalar;
-  typedef hipfftHandle  FFTW_plan;
+  typedef hipfftComplex      FFTW_scalar;
+  typedef hipfftHandle        FFTW_plan;
+
   static FFTW_plan fftw_plan_many_dft(int rank, int *n,int howmany,
-				      FFTW_scalar *in, int *inembed,
-				      int istride, int idist,
-				      FFTW_scalar *out, int *onembed,
-				      int ostride, int odist,
+				      FFTW_scalar *in, int *inembed,		
+				      int istride, int idist,		
+				      FFTW_scalar *out, int *onembed,		
+				      int ostride, int odist,		
 				      int sign, unsigned flags) {
     FFTW_plan p;
     auto rv = hipfftPlanMany(&p,rank,n,n,istride,idist,n,ostride,odist,HIPFFT_C2C,howmany);
     GRID_ASSERT(rv==HIPFFT_SUCCESS);
     return p;
-  }
+  }	  
+    
   inline static void fftw_execute_dft(const FFTW_plan p,FFTW_scalar *in,FFTW_scalar *out, int sign) {
     hipfftResult rv;
     if ( sign == forward ) rv =hipfftExecC2C(p,in,out,HIPFFT_FORWARD);
@@ -108,7 +113,9 @@ public:
     accelerator_barrier();
     GRID_ASSERT(rv==HIPFFT_SUCCESS);
   }
-  inline static void fftw_destroy_plan(const FFTW_plan p) { hipfftDestroy(p); }
+  inline static void fftw_destroy_plan(const FFTW_plan p) {
+    hipfftDestroy(p);
+  }
 };
 #endif
 
@@ -119,45 +126,53 @@ public:
   static const int backward=FFTW_BACKWARD;
   typedef cufftDoubleComplex FFTW_scalar;
   typedef cufftHandle        FFTW_plan;
+
   static FFTW_plan fftw_plan_many_dft(int rank, int *n,int howmany,
-				      FFTW_scalar *in, int *inembed,
-				      int istride, int idist,
-				      FFTW_scalar *out, int *onembed,
-				      int ostride, int odist,
+				      FFTW_scalar *in, int *inembed,		
+				      int istride, int idist,		
+				      FFTW_scalar *out, int *onembed,		
+				      int ostride, int odist,		
 				      int sign, unsigned flags) {
     FFTW_plan p;
     cufftPlanMany(&p,rank,n,n,istride,idist,n,ostride,odist,CUFFT_Z2Z,howmany);
     return p;
-  }
+  }	  
+    
   inline static void fftw_execute_dft(const FFTW_plan p,FFTW_scalar *in,FFTW_scalar *out, int sign) {
     if ( sign == forward ) cufftExecZ2Z(p,in,out,CUFFT_FORWARD);
     else                   cufftExecZ2Z(p,in,out,CUFFT_INVERSE);
     accelerator_barrier();
   }
-  inline static void fftw_destroy_plan(const FFTW_plan p) { cufftDestroy(p); }
+  inline static void fftw_destroy_plan(const FFTW_plan p) {
+    cufftDestroy(p);
+  }
 };
 template<> struct FFTW<ComplexF> {
 public:
   static const int forward=FFTW_FORWARD;
   static const int backward=FFTW_BACKWARD;
   typedef cufftComplex FFTW_scalar;
-  typedef cufftHandle  FFTW_plan;
+  typedef cufftHandle        FFTW_plan;
+
   static FFTW_plan fftw_plan_many_dft(int rank, int *n,int howmany,
-				      FFTW_scalar *in, int *inembed,
-				      int istride, int idist,
-				      FFTW_scalar *out, int *onembed,
-				      int ostride, int odist,
+				      FFTW_scalar *in, int *inembed,		
+				      int istride, int idist,		
+				      FFTW_scalar *out, int *onembed,		
+				      int ostride, int odist,		
 				      int sign, unsigned flags) {
     FFTW_plan p;
     cufftPlanMany(&p,rank,n,n,istride,idist,n,ostride,odist,CUFFT_C2C,howmany);
     return p;
-  }
+  }	  
+    
   inline static void fftw_execute_dft(const FFTW_plan p,FFTW_scalar *in,FFTW_scalar *out, int sign) {
     if ( sign == forward ) cufftExecC2C(p,in,out,CUFFT_FORWARD);
     else                   cufftExecC2C(p,in,out,CUFFT_INVERSE);
     accelerator_barrier();
   }
-  inline static void fftw_destroy_plan(const FFTW_plan p) { cufftDestroy(p); }
+  inline static void fftw_destroy_plan(const FFTW_plan p) {
+    cufftDestroy(p);
+  }
 };
 #endif
 
@@ -168,325 +183,313 @@ public:
   typedef fftw_complex FFTW_scalar;
   typedef fftw_plan    FFTW_plan;
   static FFTW_plan fftw_plan_many_dft(int rank, int *n,int howmany,
-				      FFTW_scalar *in, int *inembed,
-				      int istride, int idist,
-				      FFTW_scalar *out, int *onembed,
-				      int ostride, int odist,
+				      FFTW_scalar *in, int *inembed,		
+				      int istride, int idist,		
+				      FFTW_scalar *out, int *onembed,		
+				      int ostride, int odist,		
 				      int sign, unsigned flags) {
     return ::fftw_plan_many_dft(rank,n,howmany,in,inembed,istride,idist,out,onembed,ostride,odist,sign,flags);
-  }
+  }	  
+    
   inline static void fftw_execute_dft(const FFTW_plan p,FFTW_scalar *in,FFTW_scalar *out, int sign) {
     ::fftw_execute_dft(p,in,out);
   }
-  inline static void fftw_destroy_plan(const FFTW_plan p) { ::fftw_destroy_plan(p); }
+  inline static void fftw_destroy_plan(const FFTW_plan p) {
+    ::fftw_destroy_plan(p);
+  }
 };
 template<> struct FFTW<ComplexF> {
 public:
   typedef fftwf_complex FFTW_scalar;
   typedef fftwf_plan    FFTW_plan;
   static FFTW_plan fftw_plan_many_dft(int rank, int *n,int howmany,
-				      FFTW_scalar *in, int *inembed,
-				      int istride, int idist,
-				      FFTW_scalar *out, int *onembed,
-				      int ostride, int odist,
+				      FFTW_scalar *in, int *inembed,		
+				      int istride, int idist,		
+				      FFTW_scalar *out, int *onembed,		
+				      int ostride, int odist,		
 				      int sign, unsigned flags) {
     return ::fftwf_plan_many_dft(rank,n,howmany,in,inembed,istride,idist,out,onembed,ostride,odist,sign,flags);
-  }
+  }	  
+    
   inline static void fftw_execute_dft(const FFTW_plan p,FFTW_scalar *in,FFTW_scalar *out, int sign) {
     ::fftwf_execute_dft(p,in,out);
   }
-  inline static void fftw_destroy_plan(const FFTW_plan p) { ::fftwf_destroy_plan(p); }
+  inline static void fftw_destroy_plan(const FFTW_plan p) {
+    ::fftwf_destroy_plan(p);
+  }
 };
 #endif
 #endif
 
-struct FFTbase {
-  double        flops;
-  double        flops_call;
-  uint64_t      usec;
-  GridCartesian *_grid;
-
-  static const int forward  = FFTW_FORWARD;
-  static const int backward = FFTW_BACKWARD;
-
-  double Flops(void)  { return flops; }
-  double MFlops(void) { return flops / usec; }
-  double USec(void)   { return (double)usec; }
-
-  FFTbase(GridCartesian *grid) : _grid(grid), flops(0), flops_call(0), usec(0) {}
-};
-
-// Barrel-shift gather, FFT execute, and insert. Called by both FFT and PlannedFFT.
-// The caller is responsible for plan acquisition and destruction.
-template<class vobj>
-static void FFT_dim_execute(
-    Lattice<vobj>       &result,
-    const Lattice<vobj> &source,
-    int dim, int sign,
-    typename FFTW<typename vobj::scalar_type>::FFTW_plan p,
-    GridCartesian *grid,
-    double &flops, double &flops_call, uint64_t &usec)
-{
-  typedef typename vobj::scalar_type   scalar;
-  typedef typename vobj::scalar_object sobj;
-  typedef typename vobj::scalar_type   scalar_type;
-  typedef typename vobj::vector_type   vector_type;
-  typedef typename FFTW<scalar>::FFTW_scalar FFTW_scalar;
-
-  const int Ndim = grid->Nd();
-  int L = grid->_ldimensions[dim];
-  int G = grid->_fdimensions[dim];
-  int Ncomp = sizeof(sobj) / sizeof(scalar);
-  int64_t Nlow = 1, Nhigh = 1;
-  for (int d = 0; d < dim;      d++) Nlow  *= grid->_ldimensions[d];
-  for (int d = dim+1; d < Ndim; d++) Nhigh *= grid->_ldimensions[d];
-  int64_t Nperp = Nlow * Nhigh;
-
-  deviceVector<scalar> pgbuf(Nperp * Ncomp * G);
-  scalar *pgbuf_v = &pgbuf[0];
-  int howmany = Ncomp * Nperp;
-
-  scalar div;
-  if      (sign == FFTW_BACKWARD) div = 1.0 / G;
-  else if (sign == FFTW_FORWARD)  div = 1.0;
-  else GRID_ASSERT(0);
-
-  double t_pencil = 0, t_fft = 0, t_copy = 0, t_shift = 0;
-  double t_total  = -usecond();
-
-  result = source;
-  int pc = grid->_processor_coor[dim];
-
-  const Coordinate ldims      = grid->_ldimensions;
-  const Coordinate rdims      = grid->_rdimensions;
-  const Coordinate sdims      = grid->_simd_layout;
-  const Coordinate processors = grid->_processors;
-
-  Coordinate pgdims(Ndim);
-  pgdims[0] = G;
-  for (int d = 0, dd = 1; d < Ndim; d++)
-    if (d != dim) pgdims[dd++] = ldims[d];
-  int64_t pgvol = 1;
-  for (int d = 0; d < Ndim; d++) pgvol *= pgdims[d];
-
-  const int Nsimd = vobj::Nsimd();
-  t_pencil = -usecond();
-  for (int p_idx = 0; p_idx < processors[dim]; p_idx++) {
-    t_copy -= usecond();
-    autoView(r_v, result, AcceleratorRead);
-    accelerator_for(idx, grid->oSites(), vobj::Nsimd(), {
-#ifdef GRID_SIMT
-    {
-      int lane = acceleratorSIMTlane(Nsimd);
-#else
-    for (int lane = 0; lane < Nsimd; lane++) {
-#endif
-      Coordinate icoor, ocoor, pgcoor;
-      Lexicographic::CoorFromIndex(icoor, lane, sdims);
-      Lexicographic::CoorFromIndex(ocoor, idx,  rdims);
-      pgcoor[0] = ocoor[dim] + icoor[dim]*rdims[dim] + ((pc+p_idx)%processors[dim])*L;
-      for (int d = 0, dd = 1; d < Ndim; d++)
-        if (d != dim) { pgcoor[dd] = ocoor[d] + icoor[d]*rdims[d]; dd++; }
-      int64_t pgidx;
-      Lexicographic::IndexFromCoor(pgcoor, pgidx, pgdims);
-      vector_type *from = (vector_type *)&r_v[idx];
-      scalar_type  stmp;
-      for (int w = 0; w < Ncomp; w++) {
-        stmp = getlane(from[w], lane);
-        pgbuf_v[pgidx + w*pgvol] = stmp;
-      }
-#ifdef GRID_SIMT
-    }
-#else
-    }
-#endif
-    });
-    t_copy += usecond();
-    if (p_idx != processors[dim] - 1) {
-      Lattice<vobj> temp(grid);
-      t_shift -= usecond();
-      temp = Cshift(result, dim, L); result = temp;
-      t_shift += usecond();
-    }
-  }
-  t_pencil += usecond();
-
-  FFTW_scalar *in  = (FFTW_scalar *)pgbuf_v;
-  FFTW_scalar *out = (FFTW_scalar *)pgbuf_v;
-  t_fft = -usecond();
-  FFTW<scalar>::fftw_execute_dft(p, in, out, sign);
-  t_fft += usecond();
-
-  flops_call = 5.0 * howmany * G * log2(G);
-  usec  = t_fft;
-  flops = flops_call;
-
-  result = Zero();
-  double t_insert = -usecond();
-  {
-    autoView(r_v, result, AcceleratorWrite);
-    accelerator_for(idx, grid->oSites(), Nsimd, {
-#ifdef GRID_SIMT
-    {
-      int lane = acceleratorSIMTlane(Nsimd);
-#else
-    for (int lane = 0; lane < Nsimd; lane++) {
-#endif
-      Coordinate icoor(Ndim), ocoor(Ndim), pgcoor(Ndim);
-      Lexicographic::CoorFromIndex(icoor, lane, sdims);
-      Lexicographic::CoorFromIndex(ocoor, idx,  rdims);
-      pgcoor[0] = ocoor[dim] + icoor[dim]*rdims[dim] + pc*L;
-      for (int d = 0, dd = 1; d < Ndim; d++)
-        if (d != dim) { pgcoor[dd] = ocoor[d] + icoor[d]*rdims[d]; dd++; }
-      int64_t pgidx;
-      Lexicographic::IndexFromCoor(pgcoor, pgidx, pgdims);
-      vector_type *to  = (vector_type *)&r_v[idx];
-      scalar_type  stmp;
-      for (int w = 0; w < Ncomp; w++) {
-        stmp = pgbuf_v[pgidx + w*pgvol];
-        putlane(to[w], stmp, lane);
-      }
-#ifdef GRID_SIMT
-    }
-#else
-    }
-#endif
-    });
-  }
-  result = result * div;
-  t_insert += usecond();
-  t_total  += usecond();
-
-  std::cout << GridLogPerformance << " FFT took   " << t_total/1.0e6  << " s" << std::endl;
-  std::cout << GridLogPerformance << " FFT pencil " << t_pencil/1.0e6 << " s" << std::endl;
-  std::cout << GridLogPerformance << "  of which copy " << t_copy/1.0e6  << " s" << std::endl;
-  std::cout << GridLogPerformance << "  of which shift" << t_shift/1.0e6 << " s" << std::endl;
-  std::cout << GridLogPerformance << " FFT kernels " << t_fft/1.0e6    << " s" << std::endl;
-  std::cout << GridLogPerformance << " FFT insert  " << t_insert/1.0e6 << " s" << std::endl;
-}
-
-class FFT : public FFTbase {
+class FFT {
+private:
+    
+  double flops;
+  double flops_call;
+  uint64_t usec;
+    
 public:
-  FFT(GridCartesian *grid) : FFTbase(grid) {}
-  ~FFT() {}
+    
+  static const int forward=FFTW_FORWARD;
+  static const int backward=FFTW_BACKWARD;
+    
+  double Flops(void) {return flops;}
+  double MFlops(void) {return flops/usec;}
+  double USec(void)   {return (double)usec;}    
 
+  FFT ( GridCartesian * grid ) 
+  {
+    flops=0;
+    usec =0;
+  };
+    
+  ~FFT ( void)  {
+    //    delete sgrid;
+  }
+    
   template<class vobj>
-  void FFT_dim_mask(Lattice<vobj> &result, const Lattice<vobj> &source, Coordinate mask, int sign) {
-    const int Ndim = _grid->Nd();
+  void FFT_dim_mask(Lattice<vobj> &result,const Lattice<vobj> &source,Coordinate mask,int sign){
+
+    //    vgrid=result.Grid();
+    //    conformable(result.Grid(),vgrid);
+    //    conformable(source.Grid(),vgrid);
+    const int Ndim = source.Grid()->Nd();
     Lattice<vobj> tmp = source;
-    for (int d = 0; d < Ndim; d++) {
-      if (mask[d]) {
-        FFT_dim(result, tmp, d, sign);
-        tmp = result;
+    for(int d=0;d<Ndim;d++){
+      if( mask[d] ) {
+	FFT_dim(result,tmp,d,sign);
+	tmp=result;
       }
     }
   }
 
   template<class vobj>
-  void FFT_all_dim(Lattice<vobj> &result, const Lattice<vobj> &source, int sign) {
-    Coordinate mask(_grid->Nd(), 1);
-    FFT_dim_mask(result, source, mask, sign);
+  void FFT_all_dim(Lattice<vobj> &result,const Lattice<vobj> &source,int sign){
+    const int Ndim = source.Grid()->Nd();
+    Coordinate mask(Ndim,1);
+    FFT_dim_mask(result,source,mask,sign);
   }
 
-  template<class vobj>
-  void FFT_dim(Lattice<vobj> &result, const Lattice<vobj> &source, int dim, int sign) {
-    GRID_ASSERT(source.Grid() == _grid);
-    GRID_ASSERT(result.Grid() == _grid);
-    conformable(result.Grid(), source.Grid());
 
-    typedef typename vobj::scalar_type   scalar;
+  template<class vobj>
+  void FFT_dim(Lattice<vobj> &result,const Lattice<vobj> &source,int dim, int sign){
+    const int Ndim = source.Grid()->Nd();
+    GridBase *grid = source.Grid();
+    conformable(result.Grid(),source.Grid());
+
+    int L = grid->_ldimensions[dim];
+    int G = grid->_fdimensions[dim];
+      
+    Coordinate layout(Ndim,1);
+    
+    // Construct pencils
     typedef typename vobj::scalar_object sobj;
+    typedef typename vobj::scalar_type   scalar;
+    typedef typename vobj::scalar_type   scalar_type;
+    typedef typename vobj::vector_type   vector_type;
+      
+    //std::cout << "CPU view" << std::endl;
+    
     typedef typename FFTW<scalar>::FFTW_scalar FFTW_scalar;
     typedef typename FFTW<scalar>::FFTW_plan   FFTW_plan;
+      
+    int Ncomp = sizeof(sobj)/sizeof(scalar);
+    int64_t Nlow  = 1;
+    int64_t Nhigh = 1;
 
-    const int Ndim = _grid->Nd();
-    int G     = _grid->_fdimensions[dim];
-    int Ncomp = sizeof(sobj) / sizeof(scalar);
-    int64_t Nperp = 1;
-    for (int d = 0; d < Ndim; d++)
-      if (d != dim) Nperp *= _grid->_ldimensions[d];
-    int n[]     = {G};
-    int howmany = Ncomp * Nperp;
-
-    deviceVector<scalar> dummy(2);
-    FFTW_scalar *buf = (FFTW_scalar *)&dummy[0];
-    FFTW_plan p = FFTW<scalar>::fftw_plan_many_dft(1, n, howmany,
-                                                    buf, n, 1, G,
-                                                    buf, n, 1, G,
-                                                    sign, FFTW_ESTIMATE);
-    FFT_dim_execute(result, source, dim, sign, p, _grid, flops, flops_call, usec);
-    FFTW<scalar>::fftw_destroy_plan(p);
-  }
-};
-
-template<class vobj>
-class PlannedFFT : public FFTbase {
-private:
-  typedef typename vobj::scalar_type   scalar;
-  typedef typename vobj::scalar_object sobj;
-  typedef typename vobj::vector_type   vector_type;
-  typedef typename FFTW<scalar>::FFTW_scalar FFTW_scalar;
-  typedef typename FFTW<scalar>::FFTW_plan   FFTW_plan;
-
-  std::vector<FFTW_plan> forward_plans;
-  std::vector<FFTW_plan> backward_plans;
-
-  void PlanCreate() {
-    const int Ndim = _grid->Nd();
-    forward_plans.resize(Ndim);
-    backward_plans.resize(Ndim);
-
-    for (int d = 0; d < Ndim; d++) {
-      int G       = _grid->_fdimensions[d];
-      int Ncomp   = sizeof(sobj) / sizeof(scalar);
-      int64_t Nperp = 1;
-      for (int dd = 0; dd < Ndim; dd++)
-        if (dd != d) Nperp *= _grid->_ldimensions[dd];
-      int howmany = Ncomp * (int)Nperp;
-      int n[]     = {G};
-
-      deviceVector<scalar> dummy(2);
-      FFTW_scalar *buf = (FFTW_scalar *)&dummy[0];
-
-      forward_plans[d]  = FFTW<scalar>::fftw_plan_many_dft(1, n, howmany, buf, n, 1, G, buf, n, 1, G, FFTW_FORWARD,  FFTW_ESTIMATE);
-      backward_plans[d] = FFTW<scalar>::fftw_plan_many_dft(1, n, howmany, buf, n, 1, G, buf, n, 1, G, FFTW_BACKWARD, FFTW_ESTIMATE);
+    for(int d=0;d<dim;d++){
+      Nlow*=grid->_ldimensions[d];
     }
-  }
+    for(int d=dim+1;d<Ndim;d++){
+      Nhigh*=grid->_ldimensions[d];
+    }
+    int64_t Nperp=Nlow*Nhigh;
+    
+    deviceVector<scalar> pgbuf; // Layout is [perp][component][dim]
+    pgbuf.resize(Nperp*Ncomp*G);
+    scalar *pgbuf_v = &pgbuf[0];
+      
+    int rank = 1;  /* 1d transforms */
+    int n[] = {G}; /* 1d transforms of length G */
+    int howmany = Ncomp * Nperp;
+    int odist,idist,istride,ostride;
+    idist   = odist   = G;            /* Distance between consecutive FT's */
+    istride = ostride = 1;            /* Distance between two elements in the same FT */
+    int *inembed = n, *onembed = n;
+      
+    scalar div;
+    if ( sign == backward ) div = 1.0/G;
+    else if ( sign == forward ) div = 1.0;
+    else GRID_ASSERT(0);
 
-  void PlanDestroy() {
-    for (auto p : forward_plans)  FFTW<scalar>::fftw_destroy_plan(p);
-    for (auto p : backward_plans) FFTW<scalar>::fftw_destroy_plan(p);
-    forward_plans.clear();
-    backward_plans.clear();
-  }
+    double t_pencil=0;
+    double t_fft   =0;
+    double t_total =-usecond();
+    //    std::cout << GridLogPerformance<<"Making FFTW plan" << std::endl;
+    /*
+     *
+     */
+    FFTW_plan p;
+    {
+      FFTW_scalar *in = (FFTW_scalar *)&pgbuf_v[0];
+      FFTW_scalar *out= (FFTW_scalar *)&pgbuf_v[0];
+      p = FFTW<scalar>::fftw_plan_many_dft(rank,n,howmany,
+					   in,inembed,
+					   istride,idist,
+					   out,onembed,
+					   ostride, odist,
+					   sign,FFTW_ESTIMATE);
+    }
+      
+    // Barrel shift and collect global pencil
+    //    std::cout << GridLogPerformance<<"Making pencil" << std::endl;
+    Coordinate lcoor(Ndim), gcoor(Ndim);
+    double t_copy=0;
+    double t_shift=0;
+    t_pencil = -usecond();
+    result = source;
+    int pc = grid->_processor_coor[dim];
 
-public:
-  PlannedFFT(GridCartesian *grid) : FFTbase(grid) { PlanCreate(); }
-  ~PlannedFFT() { PlanDestroy(); }
+    const Coordinate ldims = grid->_ldimensions;
+    const Coordinate rdims = grid->_rdimensions;
+    const Coordinate sdims = grid->_simd_layout;
 
-  void FFT_dim_mask(Lattice<vobj> &result, const Lattice<vobj> &source, Coordinate mask, int sign) {
-    const int Ndim = _grid->Nd();
-    Lattice<vobj> tmp = source;
-    for (int d = 0; d < Ndim; d++) {
-      if (mask[d]) {
-        FFT_dim(result, tmp, d, sign);
-        tmp = result;
+    Coordinate processors = grid->_processors;
+    Coordinate pgdims(Ndim);
+    pgdims[0] = G;
+    for(int d=0, dd=1;d<Ndim;d++){
+      if ( d!=dim ) pgdims[dd++] = ldims[d];
+    }
+    int64_t pgvol=1;
+    for(int d=0;d<Ndim;d++) pgvol*=pgdims[d];
+    
+    const int Nsimd = vobj::Nsimd();
+    for(int p=0;p<processors[dim];p++) {
+      t_copy-=usecond();
+      autoView(r_v,result,AcceleratorRead);
+      accelerator_for(idx, grid->oSites(), vobj::Nsimd(), {
+#ifdef GRID_SIMT
+      {
+	int lane=acceleratorSIMTlane(Nsimd); // buffer lane
+#else
+      for(int lane=0;lane<Nsimd;lane++) {
+#endif
+	Coordinate icoor;
+	Coordinate ocoor;
+	Coordinate pgcoor;
+
+	Lexicographic::CoorFromIndex(icoor,lane,sdims);
+	Lexicographic::CoorFromIndex(ocoor,idx,rdims);
+
+	pgcoor[0] = ocoor[dim] + icoor[dim]*rdims[dim] + ((pc+p)%processors[dim])*L;
+	for(int d=0,dd=1;d<Ndim;d++){
+	  if ( d!=dim ) {
+	    pgcoor[dd] = ocoor[d] + icoor[d]*rdims[d];
+	    dd++;
+	  }
+	}
+
+	// Map coordinates in lattice layout to FFTW index
+	int64_t pgidx;
+	Lexicographic::IndexFromCoor(pgcoor,pgidx,pgdims);
+
+	vector_type *from = (vector_type *)&r_v[idx];
+	scalar_type stmp;
+	for(int w=0;w<Ncomp;w++){
+	  int64_t pg_idx = pgidx + w*pgvol;
+	  stmp = getlane(from[w], lane);
+	  pgbuf_v[pg_idx] = stmp;
+	}
+#ifdef GRID_SIMT
+      }
+#else
+      }
+#endif
+      });
+
+      t_copy+=usecond();
+      if (p != processors[dim] - 1) {
+	Lattice<vobj> temp(grid);
+	t_shift-=usecond();
+	temp = Cshift(result,dim,L); result = temp;
+	t_shift+=usecond();
       }
     }
-  }
+    t_pencil += usecond();
+      
+    FFTW_scalar *in = (FFTW_scalar *)pgbuf_v;
+    FFTW_scalar *out= (FFTW_scalar *)pgbuf_v;
+    t_fft = -usecond();
+    FFTW<scalar>::fftw_execute_dft(p,in,out,sign);
+    t_fft += usecond();
+    
+    // performance counting
+    flops_call = 5.0*howmany*G*log2(G);
+    usec = t_fft;
+    flops= flops_call;
 
-  void FFT_all_dim(Lattice<vobj> &result, const Lattice<vobj> &source, int sign) {
-    Coordinate mask(_grid->Nd(), 1);
-    FFT_dim_mask(result, source, mask, sign);
-  }
+    result = Zero();
+    
+    double t_insert = -usecond();
+    {
+      autoView(r_v,result,AcceleratorWrite);
+      accelerator_for(idx,grid->oSites(),Nsimd,{
+#ifdef GRID_SIMT
+      {
+	int lane=acceleratorSIMTlane(Nsimd); // buffer lane
+#else
+      for(int lane=0;lane<Nsimd;lane++) {
+#endif
+	Coordinate icoor(Ndim);
+	Coordinate ocoor(Ndim);
+	Coordinate pgcoor(Ndim);
 
-  void FFT_dim(Lattice<vobj> &result, const Lattice<vobj> &source, int dim, int sign) {
-    GRID_ASSERT(source.Grid() == _grid);
-    GRID_ASSERT(result.Grid() == _grid);
-    GRID_ASSERT((int)forward_plans.size() == _grid->Nd());
-    conformable(result.Grid(), source.Grid());
-    FFTW_plan p = (sign == forward ? forward_plans : backward_plans)[dim];
-    FFT_dim_execute(result, source, dim, sign, p, _grid, flops, flops_call, usec);
+	Lexicographic::CoorFromIndex(icoor,lane,sdims);
+	Lexicographic::CoorFromIndex(ocoor,idx,rdims);
+
+	pgcoor[0] = ocoor[dim] + icoor[dim]*rdims[dim] + pc*L;
+	for(int d=0,dd=1;d<Ndim;d++){
+	  if ( d!=dim ) {
+	    pgcoor[dd] = ocoor[d] + icoor[d]*rdims[d];
+	    dd++;
+	  }
+	}
+	// Map coordinates in lattice layout to FFTW index
+	int64_t pgidx;
+	Lexicographic::IndexFromCoor(pgcoor,pgidx,pgdims);
+
+	vector_type *to = (vector_type *)&r_v[idx];
+	scalar_type stmp;
+	for(int w=0;w<Ncomp;w++){
+	  int64_t pg_idx = pgidx + w*pgvol;
+	  stmp = pgbuf_v[pg_idx];
+	  putlane(to[w], stmp, lane);
+	}
+	
+#ifdef GRID_SIMT
+      }
+#else
+      }
+#endif
+      });
+    }
+
+    result = result*div;
+
+    t_insert +=usecond();
+    
+    // destroying plan
+    FFTW<scalar>::fftw_destroy_plan(p);
+
+    t_total +=usecond();
+
+    std::cout <<GridLogPerformance<< " FFT took   "<<t_total/1.0e6 <<" s" << std::endl;
+    std::cout <<GridLogPerformance<< " FFT pencil "<<t_pencil/1.0e6 <<" s" << std::endl;
+    std::cout <<GridLogPerformance<< "  of which copy "<<t_copy/1.0e6 <<" s" << std::endl;
+    std::cout <<GridLogPerformance<< "  of which shift"<<t_shift/1.0e6 <<" s" << std::endl;
+    std::cout <<GridLogPerformance<< " FFT kernels "<<t_fft/1.0e6 <<" s" << std::endl;
+    std::cout <<GridLogPerformance<< " FFT insert  "<<t_insert/1.0e6 <<" s" << std::endl;
+    
   }
 };
 

Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczosCoarse.h

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

Grid/lattice/Lattice_reduction.h

@@ -31,6 +31,9 @@ Author: Christoph Lehner <christoph@lhnr.de>
 #if defined(GRID_SYCL)
 #include <Grid/lattice/Lattice_reduction_sycl.h>
 #endif
+#if defined(GRID_CUDA)||defined(GRID_HIP)||defined(GRID_SYCL)
+#include <Grid/lattice/Lattice_reduction_gpu_cub.h>
+#endif
 #include <Grid/lattice/Lattice_slicesum_core.h>
 
 NAMESPACE_BEGIN(Grid);

Grid/lattice/Lattice_reduction_gpu.h

@@ -197,15 +197,12 @@ __global__ void reduceKernel(const vobj *lat, sobj *buffer, Iterator n) {
 /////////////////////////////////////////////////////////////////////////////////////////////////////////
 // Possibly promote to double and sum
 /////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-#undef GRID_REDUCTION_TIMING
-
 template <class vobj>
-inline typename vobj::scalar_objectD sumD_gpu_small(const vobj *lat, Integer osites)
+inline typename vobj::scalar_objectD sumD_gpu_small_old(const vobj *lat, Integer osites)
 {
   typedef typename vobj::scalar_objectD sobj;
   typedef decltype(lat) Iterator;
-
+  
   Integer nsimd= vobj::Nsimd();
   Integer size = osites*nsimd;
 
@@ -214,193 +211,46 @@ inline typename vobj::scalar_objectD sumD_gpu_small(const vobj *lat, Integer osi
   GRID_ASSERT(ok);
 
   Integer smemSize = numThreads * sizeof(sobj);
+  // Move out of UVM
+  // Turns out I had messed up the synchronise after move to compute stream
+  // as running this on the default stream fools the synchronise
   deviceVector<sobj> buffer(numBlocks);
   sobj *buffer_v = &buffer[0];
   sobj result;
-
-#ifdef GRID_REDUCTION_TIMING
-  RealD t_kernel = -usecond();
-#endif
   reduceKernel<<< numBlocks, numThreads, smemSize, computeStream >>>(lat, buffer_v, size);
   accelerator_barrier();
-#ifdef GRID_REDUCTION_TIMING
-  t_kernel += usecond();
-  RealD t_d2h = -usecond();
-#endif
   acceleratorCopyFromDevice(buffer_v,&result,sizeof(result));
-#ifdef GRID_REDUCTION_TIMING
-  t_d2h += usecond();
-  std::cout << GridLogDebug << "  sumD_gpu_small"
-            << " sizeof(sobj)=" << sizeof(sobj)
-            << " blocks=" << numBlocks << " threads=" << numThreads
-            << " kernel+barrier=" << t_kernel << " us"
-            << " D2H=" << t_d2h << " us" << std::endl;
-#endif
   return result;
 }
 
-// Fused pack+reduce: reads R words of each vobj at word offset 'base',
-// accumulates directly into iVector<iScalar<scalarD>,R> without staging
-// through an intermediate bundle buffer.  One HBM pass instead of three.
-template <int R, class vobj, class sobj, class Iterator>
-__device__ void packReduceBlocks(
-    const iScalar<typename vobj::vector_type> *idat,
-    sobj *g_odata, Iterator osites, int base, int words)
-{
-  constexpr Iterator nsimd = vobj::Nsimd();
-  Iterator blockSize = blockDim.x;
-
-  extern __shared__ __align__(COALESCE_GRANULARITY) unsigned char shmem_pointer[];
-  sobj *sdata = (sobj *)shmem_pointer;
-
-  Iterator tid      = threadIdx.x;
-  Iterator i        = blockIdx.x * (blockSize * 2) + threadIdx.x;
-  Iterator gridSize = blockSize * 2 * gridDim.x;
-  sobj mySum = Zero();
-
-  while (i < osites * nsimd) {
-    Iterator lane = i % nsimd;
-    Iterator ss   = i / nsimd;
-    sobj tmpD; zeroit(tmpD);
-    for (int k = 0; k < R; k++) {
-      auto w = extractLane(lane, idat[ss * words + base + k]);
-      iScalar<typename vobj::scalar_typeD> wd; wd = w;
-      tmpD._internal[k] = wd;
-    }
-    mySum += tmpD;
-
-    if (i + blockSize < osites * nsimd) {
-      lane = (i + blockSize) % nsimd;
-      ss   = (i + blockSize) / nsimd;
-      sobj tmpD2; zeroit(tmpD2);
-      for (int k = 0; k < R; k++) {
-        auto w = extractLane(lane, idat[ss * words + base + k]);
-        iScalar<typename vobj::scalar_typeD> wd; wd = w;
-        tmpD2._internal[k] = wd;
-      }
-      mySum += tmpD2;
-    }
-    i += gridSize;
-  }
-
-  reduceBlock(sdata, mySum, tid);
-  if (tid == 0) g_odata[blockIdx.x] = sdata[0];
-}
-
-template <int R, class vobj, class sobj, class Iterator>
-__global__ void packReduceKernel(
-    const iScalar<typename vobj::vector_type> *idat,
-    sobj *buffer, Iterator osites, int base, int words)
-{
-  Iterator blockSize = blockDim.x;
-
-  packReduceBlocks<R, vobj, sobj>(idat, buffer, osites, base, words);
-
-  if (gridDim.x > 1) {
-    const Iterator tid = threadIdx.x;
-    __shared__ bool amLast;
-    extern __shared__ __align__(COALESCE_GRANULARITY) unsigned char shmem_pointer[];
-    sobj *smem = (sobj *)shmem_pointer;
-
-    acceleratorFence();
-
-    if (tid == 0) {
-      unsigned int ticket = atomicInc(&retirementCount, gridDim.x);
-      amLast = (ticket == gridDim.x - 1);
-    }
-    acceleratorSynchroniseAll();
-
-    if (amLast) {
-      Iterator i = tid;
-      sobj mySum = Zero();
-      while (i < (Iterator)gridDim.x) {
-        mySum += buffer[i];
-        i += blockSize;
-      }
-      reduceBlock(smem, mySum, tid);
-      if (tid == 0) {
-        buffer[0] = smem[0];
-        retirementCount = 0;
-      }
-    }
-  }
-}
-
-template<int R, class vobj>
-inline void sumD_gpu_reduce_words(const vobj *lat, Integer osites,
-                                   typename vobj::scalar_typeD *ret_p, int base)
+template <class vobj>
+inline typename vobj::scalar_objectD sumD_gpu_large_old(const vobj *lat, Integer osites)
 {
   typedef typename vobj::vector_type  vector;
   typedef typename vobj::scalar_typeD scalarD;
-  using BundleScalarD = iVector<iScalar<scalarD>, R>;
-
-  constexpr int Nsimd = vobj::Nsimd();
-  const int words = sizeof(vobj) / sizeof(vector);
-  const iScalar<vector> *idat = (const iScalar<vector> *)lat;
-
-  Integer size = (Integer)osites * Nsimd;
-  Integer numThreads, numBlocks;
-  int ok = getNumBlocksAndThreads(size, sizeof(BundleScalarD), numThreads, numBlocks);
-  GRID_ASSERT(ok);
-
-  Integer smemSize = numThreads * sizeof(BundleScalarD);
-  deviceVector<BundleScalarD> buffer(numBlocks);
-  BundleScalarD *buffer_v = &buffer[0];
-  BundleScalarD result;
-
-#ifdef GRID_REDUCTION_TIMING
-  RealD t_kernel = -usecond();
-#endif
-  packReduceKernel<R, vobj, BundleScalarD, Integer>
-    <<<numBlocks, numThreads, smemSize, computeStream>>>
-    (idat, buffer_v, osites, base, words);
-  accelerator_barrier();
-#ifdef GRID_REDUCTION_TIMING
-  t_kernel += usecond();
-  RealD t_d2h = -usecond();
-#endif
-  acceleratorCopyFromDevice(buffer_v, &result, sizeof(result));
-#ifdef GRID_REDUCTION_TIMING
-  t_d2h += usecond();
-  std::cout << GridLogDebug << " sumD_gpu_reduce_words R=" << R
-            << " base=" << base
-            << " kernel=" << t_kernel << " D2H=" << t_d2h << " us" << std::endl;
-#endif
-
-  for (int k = 0; k < R; k++)
-    ret_p[base + k] = TensorRemove(result._internal[k]);
-}
-
-template <class vobj>
-inline typename vobj::scalar_objectD sumD_gpu_large(const vobj *lat, Integer osites)
-{
-  typedef typename vobj::vector_type    vector;
-  typedef typename vobj::scalar_typeD   scalarD;
-  typedef typename vobj::scalar_objectD sobjD;
-
-  const int words = sizeof(vobj) / sizeof(vector);
-  sobjD ret; zeroit(ret);
+  typedef typename vobj::scalar_objectD sobj;
+  sobj ret;
   scalarD *ret_p = (scalarD *)&ret;
+  
+  const int words = sizeof(vobj)/sizeof(vector);
 
-#ifdef GRID_REDUCTION_TIMING
-  RealD t_large = -usecond();
-#endif
-  int w = 0;
-  while (w + 12 <= words) { sumD_gpu_reduce_words<12>(lat, osites, ret_p, w); w += 12; }
-  while (w +  4 <= words) { sumD_gpu_reduce_words< 4>(lat, osites, ret_p, w); w +=  4; }
-  while (w      <  words) { sumD_gpu_reduce_words< 1>(lat, osites, ret_p, w); w +=  1; }
-#ifdef GRID_REDUCTION_TIMING
-  t_large += usecond();
-  std::cout << GridLogDebug << "sumD_gpu_large"
-            << " sizeof(sobjD)=" << sizeof(sobjD)
-            << " words=" << words << " total=" << t_large << " us" << std::endl;
-#endif
+  deviceVector<vector> buffer(osites);
+  vector *dat = (vector *)lat;
+  vector *buf = &buffer[0];
+  iScalar<vector> *tbuf =(iScalar<vector> *)  &buffer[0];
+  for(int w=0;w<words;w++) {
 
+    accelerator_for(ss,osites,1,{
+	buf[ss] = dat[ss*words+w];
+      });
+      
+    ret_p[w] = sumD_gpu_small_old(tbuf,osites);
+  }
   return ret;
 }
 
 template <class vobj>
-inline typename vobj::scalar_objectD sumD_gpu(const vobj *lat, Integer osites)
+inline typename vobj::scalar_objectD sumD_gpu_old(const vobj *lat, Integer osites)
 {
   typedef typename vobj::scalar_objectD sobj;
   sobj ret;
@@ -411,9 +261,9 @@ inline typename vobj::scalar_objectD sumD_gpu(const vobj *lat, Integer osites)
   int ok = getNumBlocksAndThreads(size, sizeof(sobj), numThreads, numBlocks);
   
   if ( ok ) {
-    ret = sumD_gpu_small(lat,osites);
+    ret = sumD_gpu_small_old(lat,osites);
   } else {
-    ret = sumD_gpu_large(lat,osites);
+    ret = sumD_gpu_large_old(lat,osites);
   }
   return ret;
 }
@@ -422,20 +272,20 @@ inline typename vobj::scalar_objectD sumD_gpu(const vobj *lat, Integer osites)
 // Return as same precision as input performing reduction in double precision though
 /////////////////////////////////////////////////////////////////////////////////////////////////////////
 template <class vobj>
-inline typename vobj::scalar_object sum_gpu(const vobj *lat, Integer osites) 
+inline typename vobj::scalar_object sum_gpu_old(const vobj *lat, Integer osites)
 {
   typedef typename vobj::scalar_object sobj;
   sobj result;
-  result = sumD_gpu(lat,osites);
+  result = sumD_gpu_old(lat,osites);
   return result;
 }
 
 template <class vobj>
-inline typename vobj::scalar_object sum_gpu_large(const vobj *lat, Integer osites)
+inline typename vobj::scalar_object sum_gpu_large_old(const vobj *lat, Integer osites)
 {
   typedef typename vobj::scalar_object sobj;
   sobj result;
-  result = sumD_gpu_large(lat,osites);
+  result = sumD_gpu_large_old(lat,osites);
   return result;
 }
 

Grid/lattice/Lattice_reduction_gpu_cub.h

@@ -0,0 +1,361 @@
+/*************************************************************************************
+    Grid physics library, www.github.com/paboyle/Grid
+    Source file: ./Grid/lattice/Lattice_reduction_gpu_cub.h
+    Copyright (C) 2015-2024
+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 */
+#pragma once
+
+#if defined(GRID_CUDA)
+#include <cub/cub.cuh>
+#define gpucub cub
+#define gpuError_t cudaError_t
+#define gpuSuccess cudaSuccess
+#elif defined(GRID_HIP)
+#include <hipcub/hipcub.hpp>
+#define gpucub hipcub
+#define gpuError_t hipError_t
+#define gpuSuccess hipSuccess
+#endif
+
+NAMESPACE_BEGIN(Grid);
+
+/////////////////////////////////////////////////////////////////////////////////////////////////////////
+// Unified lattice reduction using CUB (CUDA/HIP) and sycl::reduction (SYCL).
+//
+// CUDA/HIP: one accelerator_for pass per site to extract SIMD lanes and promote to sobjD,
+// then CUB/hipCUB DeviceReduce::Reduce over the resulting array.
+//
+// rocPRIM's DeviceReduce requires warpSize(64) threads per block, each holding one element
+// in shared memory: sizeof(T)*64 must fit in sharedMemPerBlock.  Large QCD objects such as
+// LatticePropagator (sobjD = 2304 bytes, 64*2304 = 147 KB) exceed this budget.
+//
+// For those types sumD_gpu_large groups the vobj's vector_type words in bundles of 4,
+// reducing each bundle as an iVector<iScalar<scalarD>,4> (64 bytes, 64*64 = 4 KB — always safe).
+// Words that do not fill a complete bundle are zero-padded.
+//
+// SYCL: sycl::reduction handles any type size through the runtime, so one path suffices.
+/////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+#if defined(GRID_CUDA) || defined(GRID_HIP)
+
+#define GRID_REDUCTION_TIMING
+
+// Direct CUB reduction on the full scalar_objectD.
+// Only safe when sizeof(sobjD)*64 <= device sharedMemPerBlock.
+// Do not call directly for large composite types (e.g. LatticePropagator).
+template<class vobj>
+inline typename vobj::scalar_objectD sumD_gpu_direct(const vobj *lat, Integer osites)
+{
+  typedef typename vobj::scalar_object  sobj;
+  typedef typename vobj::scalar_objectD sobjD;
+
+  const Integer nsimd    = vobj::Nsimd();
+  const Integer nlanes   = osites * nsimd;
+
+  deviceVector<sobjD> per_lane(nlanes);
+  sobjD *per_lane_p = &per_lane[0];
+
+#ifdef GRID_REDUCTION_TIMING
+  RealD t_for = -usecond();
+#endif
+  accelerator_for(idx, nlanes, 1, {
+    Integer ss   = idx / nsimd;
+    Integer lane = idx % nsimd;
+    sobj tmp = extractLane(lane, lat[ss]);
+    sobjD tmpD; tmpD = tmp;
+    per_lane_p[idx] = tmpD;
+  });
+#ifdef GRID_REDUCTION_TIMING
+  accelerator_barrier();
+  t_for += usecond();
+#endif
+
+  sobjD zero; zeroit(zero);
+  sobjD *d_out = static_cast<sobjD *>(acceleratorAllocDevice(sizeof(sobjD)));
+  void  *d_temp = nullptr;
+  size_t temp_bytes = 0;
+
+  gpuError_t gpuErr;
+  gpuErr = gpucub::DeviceReduce::Reduce(d_temp, temp_bytes, per_lane_p, d_out,
+                                        (int)nlanes, gpucub::Sum(), zero, computeStream);
+  if (gpuErr != gpuSuccess) {
+    std::cout << GridLogError << "sumD_gpu_direct: DeviceReduce size query failed: "
+              << gpuErr << std::endl;
+    exit(EXIT_FAILURE);
+  }
+
+  d_temp = acceleratorAllocDevice(temp_bytes);
+
+#ifdef GRID_REDUCTION_TIMING
+  RealD t_cub = -usecond();
+#endif
+  gpuErr = gpucub::DeviceReduce::Reduce(d_temp, temp_bytes, per_lane_p, d_out,
+                                        (int)nlanes, gpucub::Sum(), zero, computeStream);
+  if (gpuErr != gpuSuccess) {
+    std::cout << GridLogError << "sumD_gpu_direct: DeviceReduce failed: "
+              << gpuErr << std::endl;
+    exit(EXIT_FAILURE);
+  }
+
+  accelerator_barrier();
+#ifdef GRID_REDUCTION_TIMING
+  t_cub += usecond();
+  std::cout << GridLogMessage << "sumD_gpu_direct"
+            << " sizeof(sobjD)=" << sizeof(sobjD)
+            << " accelerator_for=" << t_for << " us"
+            << " CUB_reduce="     << t_cub  << " us" << std::endl;
+#endif
+
+  sobjD result;
+  acceleratorCopyFromDevice(d_out, &result, sizeof(sobjD));
+  acceleratorFreeDevice(d_temp);
+  acceleratorFreeDevice(d_out);
+  return result;
+}
+
+// Radix-4 word-bundle path for types too large for the direct CUB path.
+// Treats vobj as words of vector_type; groups them in bundles of 4 and reduces
+// each bundle as an iVector<iScalar<scalarD>,4> — reusing Grid's existing tensor
+// type which already has accelerator_inline operator+ and zeroit().
+// sizeof = 4 * sizeof(scalarD) <= 64 bytes; 64 * 64 = 4096 bytes, safely within
+// rocPRIM's shared-memory budget on all supported devices.
+// If words % 4 != 0, the final partial bundle is zero-padded so all unused
+// slots contribute zero to the sum.
+template<class vobj>
+inline typename vobj::scalar_objectD sumD_gpu_large(const vobj *lat, Integer osites)
+{
+  typedef typename vobj::vector_type    vector;
+  typedef typename vobj::scalar_typeD   scalarD;
+  typedef typename vobj::scalar_objectD sobjD;
+  using R4 = iVector<iScalar<scalarD>, 4>;
+
+  const int words = sizeof(vobj) / sizeof(vector);
+  const int nfull = words / 4;
+  const int rem   = words % 4;
+
+  sobjD ret; zeroit(ret);
+  scalarD *ret_p = (scalarD *)&ret;
+
+  iScalar<vector> *idat = (iScalar<vector> *)lat;
+  deviceVector<R4> buf(osites);
+  R4 *buf_p = &buf[0];
+
+  R4 zero4; zeroit(zero4);
+
+  R4 *d_out = static_cast<R4 *>(acceleratorAllocDevice(sizeof(R4)));
+  void  *d_temp = nullptr;
+  size_t temp_bytes = 0;
+
+  // Probe workspace size once — type R4 and count osites are fixed across all groups.
+  gpuError_t gpuErr;
+  gpuErr = gpucub::DeviceReduce::Reduce(d_temp, temp_bytes, buf_p, d_out,
+                                        (int)osites, gpucub::Sum(), zero4, computeStream);
+  if (gpuErr != gpuSuccess) {
+    std::cout << GridLogError << "sumD_gpu_large: DeviceReduce size query failed: "
+              << gpuErr << std::endl;
+    exit(EXIT_FAILURE);
+  }
+  d_temp = acceleratorAllocDevice(temp_bytes);
+
+#ifdef GRID_REDUCTION_TIMING
+  RealD t_for_large = 0.0, t_cub_large = 0.0;
+#endif
+
+  // Full groups of 4 words.
+  for (int g = 0; g < nfull; g++) {
+    int base = 4 * g;
+#ifdef GRID_REDUCTION_TIMING
+    t_for_large -= usecond();
+#endif
+    accelerator_for(ss, osites, 1, {
+      R4 r4;
+      r4._internal[0] = TensorRemove(Reduce(idat[ss * words + base    ]));
+      r4._internal[1] = TensorRemove(Reduce(idat[ss * words + base + 1]));
+      r4._internal[2] = TensorRemove(Reduce(idat[ss * words + base + 2]));
+      r4._internal[3] = TensorRemove(Reduce(idat[ss * words + base + 3]));
+      buf_p[ss] = r4;
+    });
+#ifdef GRID_REDUCTION_TIMING
+    accelerator_barrier();
+    t_for_large += usecond();
+    t_cub_large -= usecond();
+#endif
+    gpuErr = gpucub::DeviceReduce::Reduce(d_temp, temp_bytes, buf_p, d_out,
+                                          (int)osites, gpucub::Sum(), zero4, computeStream);
+    if (gpuErr != gpuSuccess) {
+      std::cout << GridLogError << "sumD_gpu_large: DeviceReduce failed (group "
+                << g << "): " << gpuErr << std::endl;
+      exit(EXIT_FAILURE);
+    }
+    accelerator_barrier();
+#ifdef GRID_REDUCTION_TIMING
+    t_cub_large += usecond();
+#endif
+    R4 group_result;
+    acceleratorCopyFromDevice(d_out, &group_result, sizeof(R4));
+    ret_p[base    ] = TensorRemove(group_result._internal[0]);
+    ret_p[base + 1] = TensorRemove(group_result._internal[1]);
+    ret_p[base + 2] = TensorRemove(group_result._internal[2]);
+    ret_p[base + 3] = TensorRemove(group_result._internal[3]);
+  }
+
+  // Partial last group: zero-pad unused slots so they contribute nothing to the sum.
+  if (rem > 0) {
+    int base = 4 * nfull;
+#ifdef GRID_REDUCTION_TIMING
+    t_for_large -= usecond();
+#endif
+    accelerator_for(ss, osites, 1, {
+      R4 r4; zeroit(r4);
+      for (int k = 0; k < rem; k++)
+        r4._internal[k] = TensorRemove(Reduce(idat[ss * words + base + k]));
+      buf_p[ss] = r4;
+    });
+#ifdef GRID_REDUCTION_TIMING
+    accelerator_barrier();
+    t_for_large += usecond();
+    t_cub_large -= usecond();
+#endif
+    gpuErr = gpucub::DeviceReduce::Reduce(d_temp, temp_bytes, buf_p, d_out,
+                                          (int)osites, gpucub::Sum(), zero4, computeStream);
+    if (gpuErr != gpuSuccess) {
+      std::cout << GridLogError << "sumD_gpu_large: DeviceReduce failed (partial group): "
+                << gpuErr << std::endl;
+      exit(EXIT_FAILURE);
+    }
+    accelerator_barrier();
+#ifdef GRID_REDUCTION_TIMING
+    t_cub_large += usecond();
+#endif
+    R4 partial_result;
+    acceleratorCopyFromDevice(d_out, &partial_result, sizeof(R4));
+    for (int k = 0; k < rem; k++)
+      ret_p[4 * nfull + k] = TensorRemove(partial_result._internal[k]);
+  }
+
+#ifdef GRID_REDUCTION_TIMING
+  std::cout << GridLogMessage << "sumD_gpu_large"
+            << " sizeof(sobjD)=" << sizeof(sobjD)
+            << " words=" << words << " nfull=" << nfull << " rem=" << rem
+            << " accelerator_for=" << t_for_large << " us"
+            << " CUB_reduce="     << t_cub_large  << " us" << std::endl;
+#endif
+
+  acceleratorFreeDevice(d_temp);
+  acceleratorFreeDevice(d_out);
+  return ret;
+}
+
+// Dispatch: direct CUB path for types that fit in the shared-memory budget,
+// radix-4 word-bundle path for larger types.
+// Threshold 512 bytes: 64 * 512 = 32768 bytes, within rocPRIM's
+// ROCPRIM_SHARED_MEMORY_MAX on all supported devices.
+template<class vobj>
+inline typename vobj::scalar_objectD sumD_gpu(const vobj *lat, Integer osites)
+{
+  typedef typename vobj::scalar_objectD sobjD;
+  if constexpr (sizeof(sobjD) > 512) {
+    return sumD_gpu_large(lat, osites);
+  } else {
+    return sumD_gpu_direct(lat, osites);
+  }
+}
+
+template<class vobj>
+inline typename vobj::scalar_objectD sumD_gpu_small(const vobj *lat, Integer osites)
+{
+  return sumD_gpu(lat, osites);
+}
+
+template<class vobj>
+inline typename vobj::scalar_object sum_gpu(const vobj *lat, Integer osites)
+{
+  typedef typename vobj::scalar_object sobj;
+  sobj result;
+  result = sumD_gpu(lat, osites);
+  return result;
+}
+
+template<class vobj>
+inline typename vobj::scalar_object sum_gpu_large(const vobj *lat, Integer osites)
+{
+  typedef typename vobj::scalar_object sobj;
+  sobj result;
+  result = sumD_gpu_large(lat, osites);
+  return result;
+}
+
+#endif // GRID_CUDA || GRID_HIP
+
+#if defined(GRID_SYCL)
+
+// Accumulates in sobjD throughout, fixing the precision bug in the original
+// Lattice_reduction_sycl.h which accumulated in sobj then converted at the end.
+template<class vobj>
+inline typename vobj::scalar_objectD sumD_gpu(const vobj *lat, Integer osites)
+{
+  typedef typename vobj::scalar_object  sobj;
+  typedef typename vobj::scalar_objectD sobjD;
+
+  sobjD identity; zeroit(identity);
+  sobjD ret;      zeroit(ret);
+  {
+    sycl::buffer<sobjD, 1> abuff(&ret, {1});
+    theGridAccelerator->submit([&](sycl::handler &cgh) {
+      auto Reduction = sycl::reduction(abuff, cgh, identity, std::plus<>());
+      cgh.parallel_for(sycl::range<1>{(size_t)osites},
+                       Reduction,
+                       [=](sycl::id<1> item, auto &sum) {
+                         sobj s = Reduce(lat[item[0]]);
+                         sobjD sd; sd = s;
+                         sum += sd;
+                       });
+    });
+  }
+  return ret;
+}
+
+template<class vobj>
+inline typename vobj::scalar_objectD sumD_gpu_small(const vobj *lat, Integer osites)
+{
+  return sumD_gpu(lat, osites);
+}
+
+template<class vobj>
+inline typename vobj::scalar_objectD sumD_gpu_large(const vobj *lat, Integer osites)
+{
+  return sumD_gpu(lat, osites);
+}
+
+template<class vobj>
+inline typename vobj::scalar_object sum_gpu(const vobj *lat, Integer osites)
+{
+  typedef typename vobj::scalar_object sobj;
+  sobj result;
+  result = sumD_gpu(lat, osites);
+  return result;
+}
+
+template<class vobj>
+inline typename vobj::scalar_object sum_gpu_large(const vobj *lat, Integer osites)
+{
+  return sum_gpu(lat, osites);
+}
+
+#endif // GRID_SYCL
+
+NAMESPACE_END(Grid);

Grid/lattice/Lattice_reduction_sycl.h

@@ -6,64 +6,65 @@ NAMESPACE_BEGIN(Grid);
 
 
 template <class vobj>
-inline typename vobj::scalar_objectD sumD_gpu_tensor(const vobj *lat, Integer osites)
+inline typename vobj::scalar_objectD sumD_gpu_tensor_old(const vobj *lat, Integer osites)
 {
-  typedef typename vobj::scalar_object  sobj;
+  typedef typename vobj::scalar_object sobj;
   typedef typename vobj::scalar_objectD sobjD;
 
-  sobjD identity; zeroit(identity);
-  sobjD ret;      zeroit(ret);
-  {
-    sycl::buffer<sobjD, 1> abuff(&ret, {1});
+  sobj identity; zeroit(identity);
+  sobj ret; zeroit(ret);
+  Integer nsimd= vobj::Nsimd();
+  { 
+    sycl::buffer<sobj, 1> abuff(&ret, {1});
     theGridAccelerator->submit([&](sycl::handler &cgh) {
-      auto Reduction = sycl::reduction(abuff, cgh, identity, std::plus<>());
-      cgh.parallel_for(sycl::range<1>{(size_t)osites},
-                       Reduction,
-                       [=](sycl::id<1> item, auto &sum) {
-                         sobj s = Reduce(lat[item[0]]);
-                         sobjD sd; sd = s;
-                         sum += sd;
-                       });
+      auto Reduction = sycl::reduction(abuff,cgh,identity,std::plus<>());
+      cgh.parallel_for(sycl::range<1>{osites},
+                      Reduction,
+                      [=] (sycl::id<1> item, auto &sum) {
+                        auto osite   = item[0];
+                        sum +=Reduce(lat[osite]);
+                      });
     });
   }
-  return ret;
+  sobjD dret; convertType(dret,ret);
+  return dret;
 }
 
 template <class vobj>
-inline typename vobj::scalar_objectD sumD_gpu_large(const vobj *lat, Integer osites)
+inline typename vobj::scalar_objectD sumD_gpu_large_old(const vobj *lat, Integer osites)
 {
-  return sumD_gpu_tensor(lat,osites);
+  return sumD_gpu_tensor_old(lat,osites);
 }
 template <class vobj>
-inline typename vobj::scalar_objectD sumD_gpu_small(const vobj *lat, Integer osites)
+inline typename vobj::scalar_objectD sumD_gpu_small_old(const vobj *lat, Integer osites)
 {
-  return sumD_gpu_large(lat,osites);
+  return sumD_gpu_large_old(lat,osites);
 }
 
 template <class vobj>
-inline typename vobj::scalar_objectD sumD_gpu(const vobj *lat, Integer osites)
+inline typename vobj::scalar_objectD sumD_gpu_old(const vobj *lat, Integer osites)
 {
-  return sumD_gpu_large(lat,osites);
+  return sumD_gpu_large_old(lat,osites);
 }
 
 /////////////////////////////////////////////////////////////////////////////////////////////////////////
 // Return as same precision as input performing reduction in double precision though
 /////////////////////////////////////////////////////////////////////////////////////////////////////////
 template <class vobj>
-inline typename vobj::scalar_object sum_gpu(const vobj *lat, Integer osites) 
+inline typename vobj::scalar_object sum_gpu_old(const vobj *lat, Integer osites)
 {
   typedef typename vobj::scalar_object sobj;
   sobj result;
-  result = sumD_gpu(lat,osites);
+  result = sumD_gpu_old(lat,osites);
   return result;
 }
 
 template <class vobj>
-inline typename vobj::scalar_object sum_gpu_large(const vobj *lat, Integer osites)
+inline typename vobj::scalar_object sum_gpu_large_old(const vobj *lat, Integer osites)
 {
   typedef typename vobj::scalar_object sobj;
   sobj result;
-  result = sumD_gpu_large(lat,osites);
+  result = sumD_gpu_large_old(lat,osites);
   return result;
 }
 

Grid/lattice/Lattice_slicesum_core.h

@@ -1,6 +1,7 @@
 #pragma once
 
 #if defined(GRID_CUDA)
+
 #include <cub/cub.cuh>
 #define gpucub cub
 #define gpuError_t cudaError_t
@@ -56,13 +57,8 @@ inline void sliceSumReduction_cub_small(const vobj *Data,
   //copy offsets to device
   acceleratorCopyToDeviceAsynch(&offsets[0],d_offsets,sizeof(int)*(rd+1),computeStream);
   
-#if defined(__CUDACC__) && (__CUDACC_VER_MAJOR__ >= 13)
-  #define GRID_CUB_SUM_OP ::cuda::std::plus<>{}
-#else
-  #define GRID_CUB_SUM_OP ::gpucub::Sum()
-#endif
   
-  gpuError_t gpuErr = gpucub::DeviceSegmentedReduce::Reduce(temp_storage_array, temp_storage_bytes, rb_p,d_out, rd, d_offsets, d_offsets+1, GRID_CUB_SUM_OP, zero_init, computeStream);
+  gpuError_t gpuErr = gpucub::DeviceSegmentedReduce::Reduce(temp_storage_array, temp_storage_bytes, rb_p,d_out, rd, d_offsets, d_offsets+1, ::gpucub::Sum(), zero_init, computeStream);
   if (gpuErr!=gpuSuccess) {
     std::cout << GridLogError << "Lattice_slicesum_gpu.h: Encountered error during gpucub::DeviceSegmentedReduce::Reduce (setup)! Error: " << gpuErr <<std::endl;
     exit(EXIT_FAILURE);
@@ -86,13 +82,11 @@ inline void sliceSumReduction_cub_small(const vobj *Data,
   });
   
   //issue segmented reductions in computeStream
-  gpuErr = gpucub::DeviceSegmentedReduce::Reduce(temp_storage_array, temp_storage_bytes, rb_p, d_out, rd, d_offsets, d_offsets+1, GRID_CUB_SUM_OP, zero_init, computeStream);
+  gpuErr = gpucub::DeviceSegmentedReduce::Reduce(temp_storage_array, temp_storage_bytes, rb_p, d_out, rd, d_offsets, d_offsets+1,::gpucub::Sum(), zero_init, computeStream);
   if (gpuErr!=gpuSuccess) {
     std::cout << GridLogError << "Lattice_slicesum_gpu.h: Encountered error during gpucub::DeviceSegmentedReduce::Reduce! Error: " << gpuErr <<std::endl;
     exit(EXIT_FAILURE);
   }
-
-#undef GRID_CUB_SUM_OP
   
   acceleratorCopyFromDeviceAsynch(d_out,&lvSum[0],rd*sizeof(vobj),computeStream);
   

Grid/simd/Simd.h

@@ -113,14 +113,6 @@ accelerator_inline RealD    adj(const RealD  & r){ return r; }
 accelerator_inline ComplexD adj(const ComplexD& r){ return(conjugate(r)); }
 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)); }
-#endif
-
 accelerator_inline RealF real(const RealF  & r){ return r; }
 accelerator_inline RealD real(const RealD  & r){ return r; }
 accelerator_inline RealF real(const ComplexF  & r){ return r.real(); }

Grid/threads/Accelerator.cc

@@ -3,7 +3,7 @@
 NAMESPACE_BEGIN(Grid);
 int      world_rank; // Use to control world rank for print guarding
 int      acceleratorAbortOnGpuError=1;
-uint32_t accelerator_threads=8;
+uint32_t accelerator_threads=2;
 uint32_t acceleratorThreads(void)       {return accelerator_threads;};
 void     acceleratorThreads(uint32_t t) {accelerator_threads = t;};
 

Grid/threads/Accelerator.h

@@ -96,9 +96,7 @@ void     acceleratorInit(void);
 
 #ifdef GRID_CUDA
 
-NAMESPACE_END(Grid);
 #include <cuda.h>
-NAMESPACE_BEGIN(Grid);
 
 #ifdef __CUDA_ARCH__
 #define GRID_SIMT

benchmarks/Benchmark_staggered.cc

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

benchmarks/Benchmark_staggeredF.cc

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

benchmarks/Benchmark_usqcd.cc

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

scripts/filelist

@@ -11,7 +11,7 @@ CCFILES=`find . -name '*.cc' -not -path '*/instantiation/*/*' -not -path '*/gamm
 
 ZWILS_FERMION_FILES=` find . -name '*.cc' -path '*/instantiation/*' -path '*/instantiation/ZWilsonImpl*' `
 WILS_FERMION_FILES=`  find . -name '*.cc' -path '*/instantiation/*' -path '*/instantiation/WilsonImpl*' `
-STAG_FERMION_FILES=`  find . -name '*.cc' -path '*/instantiation/*' -path '*/instantiation/StaggeredImpl*' `
+STAG_FERMION_FILES=`  find . -name '*.cc' -path '*/instantiation/*' -path '*/instantiation/Staggered*' `
 GP_FERMION_FILES=`    find . -name '*.cc' -path '*/instantiation/*' -path '*/instantiation/Gparity*' `
 ADJ_FERMION_FILES=`   find . -name '*.cc' -path '*/instantiation/*' -path '*/instantiation/WilsonAdj*' `
 TWOIND_FERMION_FILES=`find . -name '*.cc' -path '*/instantiation/*' -path '*/instantiation/WilsonTwoIndex*'`

skills/gpu-memory-performance.md

@@ -1,205 +0,0 @@
----
-name: gpu-memory-performance
-description: Diagnose and fix GPU memory bandwidth and occupancy problems in Grid HPC kernels — acceleratorThreads() pitfalls, LambdaApply thread mapping, coalescedRead/Write idiom, when to use accelerator_for vs a hand-rolled __global__ kernel, and fused vs staged HBM access patterns.
-user-invocable: true
-allowed-tools:
-  - Read
-  - Bash(grep -r)
----
-
-# GPU Memory Performance in Grid
-
-## Nsimd on GPU builds
-
-With `GEN_SIMD_WIDTH=64B` (the typical production setting), `Nsimd` is **not 1**:
-
-| Scalar type | `sizeof` | `Nsimd = 64 / sizeof` |
-|---|---|---|
-| `ComplexD` | 16 B | **4** |
-| `ComplexF` | 8 B  | **8** |
-| `RealD`    | 8 B  | **8** |
-| `RealF`    | 4 B  | **16** |
-
-So for `LatticePropagatorD` (scalar type `ComplexD`), `Nsimd=4` and the SIMD lane runs `threadIdx.x ∈ {0,1,2,3}`.  True `Nsimd=1` scalar-GPU builds are the exception, not the rule.
-
-## The acceleratorThreads() Trap
-
-`acceleratorThreads()` is a runtime-settable global (default **8**) that controls the `blockDim.y` of every `accelerator_for` launch.  It is NOT the SIMD width — it is the number of sites processed per block in the y-dimension.
-
-```cpp
-// Grid/threads/Accelerator.cc
-uint32_t accelerator_threads = 8;
-```
-
-With `accelerator_for(ss, osites, nsimd, ...)`, the launch is:
-
-```
-dim3 threads(nsimd, acceleratorThreads(), 1)
-dim3 blocks ((osites + acceleratorThreads() - 1) / acceleratorThreads(), 1, 1)
-```
-
-Total threads per block = `Nsimd × acceleratorThreads()`.  With `GEN_SIMD_WIDTH=64B` and `Nsimd=8` (fp32 / ComplexF):
-
-| `acceleratorThreads()` | threads/block | AMD wavefront | note |
-|---|---|---|---|
-| 2 (original default) | 16 | 25% | sub-wavefront, poor |
-| 4 | 32 | 50% | half wavefront |
-| **8 (current default)** | **64** | **100%** | **one full wavefront — Dslash sweet spot** |
-| 16 | 128 | 200% | two wavefronts/block; register-pressure cliff for heavy kernels |
-| 32 | 256 | 400% | severe register spill for stencil kernels |
-
-**Why 8 and not higher?**  Compute-heavy kernels like the Domain Wall Dslash carry many live registers per thread (spinors + gauge links + projections).  Doubling `acceleratorThreads` from 8→16 doubles the register demand per block, which on AMD GFX90A triggers a hard occupancy cliff: `Benchmark_dwf_fp32` drops from 1.7 TF/s (nt=8) to ~300 GF/s (nt=16).  The sweet spot is one full wavefront per block, which with `Nsimd=8` (fp32) means `nt=8`.
-
-For fp64 work (`Nsimd=4`), `nt=8` gives 32 threads = half a wavefront (AMD pads to 64 with idle lanes).  Kernels that are not register-limited (e.g. simple lattice arithmetic) can benefit from `--accelerator-threads 16` at runtime.  Reduction kernels bypass `acceleratorThreads()` entirely via `getNumBlocksAndThreads`.
-
-**Why was the default ever 2?**  Before the `threadIdx.x`/`threadIdx.y` remap (see LambdaApply section below), the site index lived in `threadIdx.x` — the fast, coalescing dimension.  Increasing `acceleratorThreads` widened the block in the *site* direction, so adjacent threads in a warp hit adjacent sites, each stride `sizeof(vobj)` apart in AoS memory — breaking coalescing.  On early NVIDIA ports with `Nsimd≈8`, `nt=2` gave 16 threads = 50% of a 32-thread warp; NVIDIA recovers this via multiple concurrent blocks per SM, so occupancy was barely tolerable.  AMD has no such multiplier when blocks are already sub-wavefront.  After the remap put the site index in `threadIdx.y` and the SIMD lane in `threadIdx.x`, coalescing became independent of `acceleratorThreads`, removing the constraint.
-
-**Diagnostic**: observed bandwidth << peak, kernel time >> expected from data volume.  Check with `--accelerator-threads 32` at runtime.  A large speedup confirms occupancy starvation.
-
-**Fix options** (in order of preference):
-1. Kernel needs its own thread count — use `getNumBlocksAndThreads` and launch a `__global__` kernel directly (see below).
-2. Temporarily acceptable: set `--accelerator-threads 32` at the application level.  Note this affects every `accelerator_for` site in the binary.
-
-## LambdaApply Thread Mapping
-
-`accelerator_for` and `accelerator_for2d` go through `LambdaApply`:
-
-```cpp
-// HIP/CUDA LambdaApply kernel:
-uint64_t x = threadIdx.y + blockDim.y * blockIdx.x;  // iter1 (site index)
-uint64_t y = threadIdx.z + blockDim.z * blockIdx.y;  // iter2
-uint64_t z = threadIdx.x;                             // lane (SIMD lane)
-Lambda(x, y, z);
-```
-
-`threadIdx.x` is the **fast** (lane) dimension — consecutive thread IDs within a warp/wavefront correspond to consecutive lane values on the **same** site, not consecutive sites.
-
-With `GEN_SIMD_WIDTH=64B` and `Nsimd=4` (PropagatorD), a 64-thread AMD wavefront contains 64/4 = 16 sites, each processed by 4 lanes.  Adjacent threads within the wavefront read different lanes of the same site — this is a broadcast pattern (hardware handles this efficiently), not a stride.  The stride between consecutive *sites* (`sizeof(vobj)`) only appears between groups of `Nsimd` threads, spaced `Nsimd` apart in threadIdx.x — not between adjacent threads within a warp.
-
-## coalescedRead / coalescedWrite
-
-These are Grid's canonical way to read/write one SIMD lane from a vector type inside a `GRID_SIMT` kernel:
-
-```cpp
-// accelerator_for(ss, osites, Nsimd, {
-//   lane = acceleratorSIMTlane(Nsimd) = threadIdx.x  ∈ {0..Nsimd-1}
-auto scalar_val = coalescedRead(field[ss]);       // extractLane(lane, field[ss])
-coalescedWrite(field[ss], scalar_val);             // insertLane(lane, field[ss], scalar_val)
-```
-
-For `vobj` aggregate types, `coalescedRead` calls `extractLane(lane, vobj)` which recurses through the tensor hierarchy and returns `vobj::scalar_object`.
-
-For `vsimd` (raw SIMD vector) types, it casts to `scalar_type*` and indexes with `lane`.
-
-## Coalescing the Iteration Structure
-
-For an AoS input array where each site is `words` 16-byte elements, adjacent threads reading the same site's consecutive words achieve coalesced access:
-
-```cpp
-// Good: k varies across threads in a block → consecutive 16-byte reads
-accelerator_for2d(k, R, ss, osites, Nsimd, {
-    coalescedWrite(out[ss]._internal[k],
-                   coalescedRead(idat[ss * words + base + k]));
-});
-// dim3(Nsimd, nt, 1): threadIdx.y = k (consecutive words, coalesced)
-//                      threadIdx.x = lane (SIMD sub-lane, coalesced for Nsimd>1)
-```
-
-```cpp
-// Bad: each thread reads all R words of its site serially
-accelerator_for(ss, osites, 1, {
-    Bundle b;
-    for (int k = 0; k < R; k++)
-        b._internal[k] = idat[ss * words + base + k];  // serial, not coalesced across threads
-    out[ss] = b;                                         // bulk struct write
-});
-```
-
-The bad pattern also accumulates a large struct in registers (192 bytes for R=12), increasing register pressure and reducing occupancy further.
-
-## When to Use a __global__ Kernel Instead of accelerator_for
-
-`accelerator_for` is correct for site-parallel work where `acceleratorThreads()` is tuned appropriately.  Use a direct `__global__` kernel when:
-
-- The kernel requires a **specific thread count** for correctness or performance (reductions, shared-memory algorithms).
-- The optimal thread count depends on `sizeof(sobj)` and `sharedMemPerBlock`, not on a runtime global.
-- You need the retirement-count pattern for cross-block final reduction.
-
-Pattern: use `getNumBlocksAndThreads` to pick `numThreads` and `numBlocks`:
-
-```cpp
-Integer numThreads, numBlocks;
-int ok = getNumBlocksAndThreads(n, sizeof(sobj), numThreads, numBlocks);
-// starts at warpSize (32/64), doubles while 2*threads*sizeof(sobj) < sharedMemPerBlock
-// gives 64–256 threads/block → correct occupancy independent of acceleratorThreads()
-Integer smemSize = numThreads * sizeof(sobj);
-myKernel<<<numBlocks, numThreads, smemSize, computeStream>>>(args...);
-```
-
-Grid's `reduceKernel` uses this pattern and achieves ~400 GB/s on MI250X.
-
-## Fused vs Staged HBM Access
-
-A staged pack+reduce reads the data **three times**:
-
-```
-pack kernel:   reads vobj array (N bytes), writes bundle buffer (N bytes)
-reduce kernel: reads bundle buffer (N bytes), writes tiny result buffer
-```
-
-Total HBM: 3N bytes for N bytes of useful input.
-
-A fused kernel reads the data **once**:
-
-```
-packReduceKernel: reads R words of vobj array (N bytes), reduces in-place
-```
-
-Total HBM: N bytes.  Register pressure increases (R words held per thread) but the 3× HBM saving dominates for large objects.
-
-The fused pattern in Grid's `sumD_gpu_reduce_words<R>`:
-
-```cpp
-template <int R, class vobj, class sobj, class Iterator>
-__device__ void packReduceBlocks(
-    const iScalar<typename vobj::vector_type> *idat,
-    sobj *g_odata, Iterator osites, int base, int words)
-{
-  // sobj = iVector<iScalar<scalarD>, R>  (R double-precision scalars per site)
-  constexpr Iterator nsimd = vobj::Nsimd();
-  ...
-  while (i < osites * nsimd) {
-    Iterator lane = i % nsimd;
-    Iterator ss   = i / nsimd;
-    sobj tmpD; zeroit(tmpD);
-    for (int k = 0; k < R; k++) {
-      auto w = extractLane(lane, idat[ss * words + base + k]);
-      iScalar<typename vobj::scalar_typeD> wd; wd = w;   // float→double promotion
-      tmpD._internal[k] = wd;
-    }
-    mySum += tmpD;
-    ...
-  }
-  reduceBlock(sdata, mySum, tid);
-}
-```
-
-Launched with `getNumBlocksAndThreads` → 128 threads/block for R=12 (`BundleScalarD`=192 B, sharedMem=64 KB) → correct occupancy without depending on `acceleratorThreads()`.
-
-## Observed Numbers on MI250X (32^4 LatticePropagatorD, Nsimd=4, GEN_SIMD_WIDTH=64B)
-
-| Configuration | pack µs/group | reduce µs/group | total µs | GB/s |
-|---|---|---|---|---|
-| acceleratorThreads=2 (8 threads/block), staged | 10,080 | 470 | 126,909 | 50 |
-| acceleratorThreads=16 (64 threads/block), staged | 342 | 310 | 8,251 | 297 |
-| acceleratorThreads=16, fused (128 threads/block via getNumBlocksAndThreads) | — | 349 | 4,584 | 546 |
-
-The fused kernel at 349 µs/group reads 201 MB at 576 GB/s — 36% of MI250X HBM peak.  The remaining gap from peak is the in-kernel serial loop over R=12 words and the 12 serial kernel launches.
-
-## Quick Checklist When a Kernel Is Slow
-
-1. **Check Nsimd**: `GEN_SIMD_WIDTH=64B` → Nsimd=4 (ComplexD), 8 (ComplexF).  Total threads/block = `Nsimd × acceleratorThreads()`.  With old default nt=2 and Nsimd=4: 8 threads = 12.5% of AMD wavefront.
-2. Check threads per block: for `accelerator_for` kernels use `--accelerator-threads 32` and measure; a large speedup confirms occupancy starvation.
-3. Check for bulk struct accumulation in registers (`Bundle b; for(...) b._internal[k] = ...;`).  Replace with per-element writes via `coalescedWrite`.
-4. Check for staged HBM access (pack → buffer → reduce).  Count the passes; fuse if ≥ 2 passes over the same data.
-5. For reduction kernels, always use `getNumBlocksAndThreads` rather than `accelerator_for` so thread count is independent of `acceleratorThreads()`.

skills/mpi-heterogeneous.md

@@ -1,6 +1,6 @@
 ---
 name: mpi-heterogeneous
-description: Diagnose and work around MPI correctness bugs on heterogeneous (CPU+GPU) systems — device buffer aliasing in MPI_Sendrecv, AARCH64 PLT corruption from libfabric, topology-dependent allreduce hangs, mixed-ABI HIP runtime from wrong GTL library (Frontier/ROCm), and deterministic point-to-point reduction trees as a replacement for MPI_Allreduce.
+description: Diagnose and work around MPI correctness bugs on heterogeneous (CPU+GPU) systems — device buffer aliasing in MPI_Sendrecv, AARCH64 PLT corruption from libfabric, topology-dependent allreduce hangs, and deterministic point-to-point reduction trees as a replacement for MPI_Allreduce.
 user-invocable: true
 allowed-tools:
   - Read
@@ -110,51 +110,6 @@ void GlobalSumP2P(double *data, int count, MPI_Comm comm) {
 
 Grid reference: `USE_GRID_REDUCTION` macro in `Grid/communicator/Communicator_mpi3.cc`.
 
-## Bug Class 4: Mixed HIP ABI from Wrong GTL Library (Frontier / ROCm)
-
-**Symptom**: `HIPFFT_PARSE_ERROR` (error code 12) returned by `hipfftPlanMany` / `hipfftMakePlanMany` / `hipfftPlan1d` for FFT sizes G < 32, but G ≥ 32 succeeds. The failure only occurs with an empty rocFFT kernel cache (`~/.cache/rocfft`); a warm cache may mask it. Host-side operations and GPU kernels that do not invoke rocFFT JIT work correctly.
-
-**Root cause — mixed HIP ABI**: rocFFT uses JIT compilation (via `libamd_comgr`) for small transforms (G < 32); for G ≥ 32 it uses pre-compiled device code bundled in the library, so the JIT path is never exercised. When two HIP runtime versions are loaded in the same process — e.g. `libamdhip64.so.7` (ROCm 7) and `libamdhip64.so.6` (ROCm 6) — the rocFFT JIT cannot complete successfully.
-
-The hidden source of the old library is the Cray MPI GPU Transport Layer. On Frontier, `cray-mpich`'s `libmpi_gtl_hsa.so` may be compiled against `libamdhip64.so.6` (ROCm 6 ABI) even when the loaded ROCm module is 7.0.2. Because `LD_LIBRARY_PATH` picks up the GTL directory before the ROCm 7 library directory, `libamdhip64.so.6` is pulled in first, and both ABI versions end up resident in the process.
-
-**Diagnosis**:
-```bash
-# Check which libamdhip64 versions are actually linked into your binary at runtime
-ldd --verbose ./your_binary 2>&1 | grep amdhip
-# Bad output — two different .so versions:
-#   libamdhip64.so.6 => /opt/rocm-6.4.2/lib/libamdhip64.so.6
-#   libamdhip64.so.7 => /opt/rocm-7.0.2/lib/libamdhip64.so.7
-# Good output — only one:
-#   libamdhip64.so.7 => /opt/rocm-7.0.2/lib/libamdhip64.so.7
-```
-
-If two versions appear, the problem is the GTL/LD_LIBRARY_PATH ordering.
-
-**Fix — correct module stack and LD_LIBRARY_PATH ordering (Frontier)**:
-```bash
-module load cce/21.0.0
-module load cpe/26.03
-module load rocm/7.0.2
-# Prepend CRAY_LD_LIBRARY_PATH so the ROCm-7-aware GTL is found first
-export LD_LIBRARY_PATH=$CRAY_LD_LIBRARY_PATH:$LD_LIBRARY_PATH
-# Ensure ROCm 7 LLVM libs (needed by libamd_comgr JIT) are on the path
-export LD_LIBRARY_PATH=/opt/rocm-7.0.2/lib/llvm/lib/:$LD_LIBRARY_PATH
-```
-
-The critical step is prepending `CRAY_LD_LIBRARY_PATH`: this ensures the GTL library built against the ROCm 7 ABI is resolved before any older version that may appear further down `LD_LIBRARY_PATH`. Without this step, a stale symlink or directory ordering can silently load the wrong `libmpi_gtl_hsa.so`.
-
-**Reproducer**: `tests/debug/Test_hipfft_repro.cc` — standalone hipFFT test (no Grid headers) that sweeps G and howmany values matching realistic Grid lattice geometries. Compile with:
-```bash
-hipcc -o Test_hipfft_repro Test_hipfft_repro.cc -lhipfft
-rm -rf ~/.cache/rocfft   # empty cache required to trigger JIT path
-./Test_hipfft_repro
-```
-
-**Reference**: `systems/WorkArounds.txt`, Frontier section — GPU mapping, XPMEM, and `FI_MR_CACHE_MONITOR=disabled` settings for Frontier are documented there.
-
-**Systems affected**: Frontier (ORNL, MI250X). Likely applies to any Cray PE system where the loaded `cray-mpich` GTL was compiled against an older ROCm ABI than the runtime ROCm module. LumiG (CSC, MI250X) uses the same Cray PE and may exhibit the same issue.
-
 ## Compile-Time Guard Structure
 
 Recommended macro structure to switch between the workaround paths:

systems/Frontier/benchmarks/Benchmark_usqcd.csv

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

systems/Frontier/config-command

@@ -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,13 +9,12 @@
 --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 " \
- LDFLAGS="-L${ROCM_PATH}/lib -L${MPICH_DIR}/lib -lmpi -lmpi_gtl_hsa -lhipblas -lrocblas -lhipfft -lamdhip64"
+CXXFLAGS="-fPIC -I${ROCM_PATH}/include/ -I${MPICH_DIR}/include -L/lib64 " \
+ LDFLAGS="-L/lib64 -L${ROCM_PATH}/lib -L${MPICH_DIR}/lib -lmpi -L${CRAY_MPICH_ROOTDIR}/gtl/lib -lmpi_gtl_hsa -lhipblas -lrocblas -lhipfft"
 
 
 

systems/Frontier/sourceme.sh

@@ -1,14 +1,28 @@
 
 echo spack
-. /autofs/nccs-svm1_home1/paboyle/spack/share/spack/setup-env.sh
+. /autofs/nccs-svm1_home1/paboyle/Crusher/Grid/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/15.0.1
+module load amd/7.0.2
+#module load amd/7.1.1
+#module load rocm/7.2.0
+#module load rocm/6.4.2
+module load cray-fftw
+module load craype-accel-amd-gfx90a
 
-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
+#Ugly hacks to get down level software working on current system
+export LD_LIBRARY_PATH=/opt/cray/libfabric/1.20.1/lib64/:$LD_LIBRARY_PATH
+export LD_LIBRARY_PATH=/opt/gcc/mpfr/3.1.4/lib:$LD_LIBRARY_PATH
+export LD_LIBRARY_PATH=`pwd`/:$LD_LIBRARY_PATH
+export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$HOME/LD_PATH/
+
+#echo spack load c-lime
+#spack load c-lime
+#module load emacs 
+##module load PrgEnv-gnu
+##module load cray-mpich
+##module load cray-fftw
+##module load craype-accel-amd-gfx90a
+##export LD_LIBRARY_PATH=/opt/gcc/mpfr/3.1.4/lib:$LD_LIBRARY_PATH
+#Hack for lib
+##export LD_LIBRARY_PATH=`pwd`/:$LD_LIBRARY_PATH

systems/Perlmutter/config-command

@@ -1,12 +1,12 @@
 DIR=`pwd`
-
+PREFIX=$HOME/DDHMC/Grid/systems/Prerequisites/install/
 ../../configure \
     --enable-comms=mpi \
     --enable-simd=GPU \
     --enable-shm=nvlink \
     --enable-gen-simd-width=64 \
-    --with-gmp=$GMP \
-    --with-mpfr=$MPFR \
+    --with-gmp=$PREFIX \
+    --with-mpfr=$PREFIX \
     --enable-accelerator=cuda \
     --disable-fermion-reps \
     --disable-unified \

systems/Perlmutter/sourceme.sh

@@ -1,6 +1,4 @@
-export CRAY_ACCEL_TARGET=nvidia80
-source /global/homes/p/pboyle/spack/share/spack/setup-env.sh
-export MPFR=`spack find --paths mpfr | grep mpfr | cut -c 13-`
-export GMP=`spack find --paths gmp   | grep gmp  | cut -c 12-`
 
-module load PrgEnv-gnu cpe-cuda cudatoolkit/12.0
+export CRAY_ACCEL_TARGET=nvidia80
+
+module load PrgEnv-gnu cpe-cuda cudatoolkit/11.4

systems/mac-arm/config-command-homebrew

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

systems/mac-arm/config-command-mpi

@@ -3,10 +3,7 @@
 CXX=mpicxx ../../configure \
 	   --enable-simd=GEN \
 	   --enable-comms=mpi-auto \
-	   --enable-Sp=no \
-	   --disable-fermion-reps \
-	   --disable-gparity \
-	   --with-fftw=$FFTW \
+	   --enable-Sp=yes \
 	   --enable-unified=yes \
 	   --prefix /Users/peterboyle/QCD/vtk/Grid/install \
 	   --with-lime=$CLIME \

systems/mac-arm/sourceme-homebrew.sh

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

systems/mac-arm/sourceme.sh

@@ -1,11 +0,0 @@
-source /Users/peterboyle/QCD//Spack/spack//share/spack/setup-env.sh
-
-export FFTW=`spack find --paths fftw    | grep ^fftw   | awk '{print $2}' `
-#export HDF5=`spack find --paths hdf5+cxx   | grep ^hdf5   | awk '{print $2}' `
-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}' `
-
-export LD_LIBRARY_PATH=$MPFR/lib:$LD_LIBRARY_PATH
-export LD_LIBRARY_PATH=$GMP/lib:$LD_LIBRARY_PATH

tests/core/Test_planned_fft.cc

@@ -1,321 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid
-
-    Source file: ./tests/core/Test_planned_fft.cc
-
-    Copyright (C) 2015
-
-Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
-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 */
-#include <Grid/Grid.h>
-
-using namespace Grid;
-
-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;
-
-  Coordinate latt_size   = GridDefaultLatt();
-  Coordinate simd_layout = GridDefaultSimd(Nd,vComplexD::Nsimd());
-  Coordinate mpi_layout  = GridDefaultMpi();
-
-  int vol = 1;
-  for(int d=0;d<latt_size.size();d++) vol *= latt_size[d];
-
-  GridCartesian         GRID(latt_size,simd_layout,mpi_layout);
-  GridRedBlackCartesian RBGRID(&GRID);
-
-  LatticeComplexD     one(&GRID);
-  LatticeComplexD      zz(&GRID);
-  LatticeComplexD       C(&GRID);
-  LatticeComplexD  Ctilde(&GRID);
-  LatticeComplexD  Cref  (&GRID);
-  LatticeComplexD  Csav  (&GRID);
-  LatticeComplexD    coor(&GRID);
-
-  LatticeSpinMatrixD    S(&GRID);
-  LatticeSpinMatrixD    Stilde(&GRID);
-
-  Coordinate p({1,3,2,3});
-
-  one = ComplexD(1.0,0.0);
-  zz  = ComplexD(0.0,0.0);
-  ComplexD ci(0.0,1.0);
-
-  std::cout<<"*************************************************"<<std::endl;
-  std::cout<<"Testing Fourier form of known plane wave         "<<std::endl;
-  std::cout<<"*************************************************"<<std::endl;
-  C=Zero();
-  for(int mu=0;mu<4;mu++){
-    RealD TwoPiL =  M_PI * 2.0/ latt_size[mu];
-    LatticeCoordinate(coor,mu);
-    C = C + (TwoPiL * p[mu]) * coor;
-  }
-  C = exp(C*ci);
-  Csav = C;
-  S=Zero();
-  S = S+C;
-
-  // PlannedFFT is templated on the lattice element type (vector_object), not the Lattice<> itself.
-  PlannedFFT<LatticeComplexD::vector_object>    theFFT(&GRID);
-  PlannedFFT<LatticeSpinMatrixD::vector_object> theFFT_spin(&GRID);
-
-  Ctilde=C;
-  std::cout<<" Benchmarking PlannedFFT of LatticeComplex  "<<std::endl;
-  theFFT.FFT_dim(Ctilde,Ctilde,0,FFTbase::forward); std::cout << theFFT.MFlops()<<" Mflops "<<std::endl;
-  theFFT.FFT_dim(Ctilde,Ctilde,1,FFTbase::forward); std::cout << theFFT.MFlops()<<" Mflops "<<std::endl;
-  theFFT.FFT_dim(Ctilde,Ctilde,2,FFTbase::forward); std::cout << theFFT.MFlops()<<" Mflops "<<std::endl;
-  theFFT.FFT_dim(Ctilde,Ctilde,3,FFTbase::forward); std::cout << theFFT.MFlops()<<" Mflops "<<std::endl;
-
-  TComplexD cVol;
-  cVol()()() = vol;
-
-  Cref=Zero();
-  pokeSite(cVol,Cref,p);
-
-  Cref=Cref-Ctilde;
-  std::cout << "diff scalar "<<norm2(Cref) << std::endl;
-
-  C=Csav;
-  theFFT.FFT_all_dim(Ctilde,C,FFTbase::forward);
-  theFFT.FFT_all_dim(Cref,Ctilde,FFTbase::backward);
-
-  std::cout << norm2(C) << " " << norm2(Ctilde) << " " << norm2(Cref)<< " vol " << vol<< std::endl;
-
-  Cref= Cref - C;
-  std::cout << " invertible check " << norm2(Cref)<<std::endl;
-
-  Stilde=S;
-  std::cout<<" Benchmarking PlannedFFT of LatticeSpinMatrix  "<<std::endl;
-  theFFT_spin.FFT_dim(Stilde,Stilde,0,FFTbase::forward); std::cout << theFFT_spin.MFlops()<<" mflops "<<std::endl;
-  theFFT_spin.FFT_dim(Stilde,Stilde,1,FFTbase::forward); std::cout << theFFT_spin.MFlops()<<" mflops "<<std::endl;
-  theFFT_spin.FFT_dim(Stilde,Stilde,2,FFTbase::forward); std::cout << theFFT_spin.MFlops()<<" mflops "<<std::endl;
-  theFFT_spin.FFT_dim(Stilde,Stilde,3,FFTbase::forward); std::cout << theFFT_spin.MFlops()<<" mflops "<<std::endl;
-
-  SpinMatrixD Sp;
-  Sp = Zero(); Sp = Sp+cVol;
-
-  S=Zero();
-  pokeSite(Sp,S,p);
-
-  S= S-Stilde;
-  std::cout << "diff FT[SpinMat] "<<norm2(S) << std::endl;
-
-  std::vector<int> seeds({1,2,3,4});
-  GridSerialRNG          sRNG;  sRNG.SeedFixedIntegers(seeds);
-  GridParallelRNG          pRNG(&GRID);
-  pRNG.SeedFixedIntegers(seeds);
-
-  LatticeGaugeFieldD Umu(&GRID);
-  SU<Nc>::ColdConfiguration(pRNG,Umu);
-
-  ////////////////////////////////////////////////////
-  // Wilson test
-  ////////////////////////////////////////////////////
-  {
-    LatticeFermionD    src(&GRID); gaussian(pRNG,src);
-    LatticeFermionD    tmp(&GRID);
-    LatticeFermionD    ref(&GRID);
-
-    RealD mass=0.01;
-    WilsonFermionD Dw(Umu,GRID,RBGRID,mass);
-
-    Dw.M(src,tmp);
-
-    std::cout << "Dw src = " <<norm2(src)<<std::endl;
-    std::cout << "Dw tmp = " <<norm2(tmp)<<std::endl;
-
-    Dw.FreePropagator(tmp,ref,mass);
-
-    std::cout << "Dw ref = " <<norm2(ref)<<std::endl;
-
-    ref = ref - src;
-    std::cout << "Dw ref-src = " <<norm2(ref)<<std::endl;
-  }
-
-  ////////////////////////////////////////////////////
-  // Dwf matrix — verify Fourier representation using PlannedFFT<LatticeFermionD>
-  ////////////////////////////////////////////////////
-  {
-    std::cout<<"****************************************"<<std::endl;
-    std::cout<<"Testing Fourier representation of Ddwf"<<std::endl;
-    std::cout<<"****************************************"<<std::endl;
-
-    const int Ls=16;
-    const int sdir=0;
-    RealD mass=0.01;
-    RealD M5  =1.0;
-    Gamma G5(Gamma::Algebra::Gamma5);
-
-    GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,&GRID);
-    GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,&GRID);
-
-    DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,GRID,RBGRID,mass,M5);
-
-    GridParallelRNG          RNG5(FGrid);      RNG5.SeedFixedIntegers(seeds);
-    LatticeFermionD    src5(FGrid); gaussian(RNG5,src5);
-    LatticeFermionD    src5_p(FGrid);
-    LatticeFermionD    result5(FGrid);
-    LatticeFermionD    ref5(FGrid);
-    LatticeFermionD    tmp5(FGrid);
-
-    Ddwf.M(src5,tmp5);
-    ref5 = tmp5;
-
-    PlannedFFT<LatticeFermionD::vector_object> theFFT5(FGrid);
-
-    theFFT5.FFT_dim(result5,tmp5,1,FFTbase::forward); tmp5 = result5;
-    std::cout<<"Fourier xformed Ddwf 1 "<<norm2(result5)<<std::endl;
-    theFFT5.FFT_dim(result5,tmp5,2,FFTbase::forward); tmp5 = result5;
-    std::cout<<"Fourier xformed Ddwf 2 "<<norm2(result5)<<std::endl;
-    theFFT5.FFT_dim(result5,tmp5,3,FFTbase::forward); tmp5 = result5;
-    std::cout<<"Fourier xformed Ddwf 3 "<<norm2(result5)<<std::endl;
-    theFFT5.FFT_dim(result5,tmp5,4,FFTbase::forward);
-    std::cout<<"Fourier xformed Ddwf 4 "<<norm2(result5)<<std::endl;
-    result5 = result5*ComplexD(::sqrt(1.0/vol),0.0);
-
-    std::cout<<"Fourier xformed Ddwf "<<norm2(result5)<<std::endl;
-
-    tmp5 = src5;
-    theFFT5.FFT_dim(src5_p,tmp5,1,FFTbase::forward); tmp5 = src5_p;
-    theFFT5.FFT_dim(src5_p,tmp5,2,FFTbase::forward); tmp5 = src5_p;
-    theFFT5.FFT_dim(src5_p,tmp5,3,FFTbase::forward); tmp5 = src5_p;
-    theFFT5.FFT_dim(src5_p,tmp5,4,FFTbase::forward); src5_p = src5_p*ComplexD(::sqrt(1.0/vol),0.0);
-
-    std::cout<<"Fourier xformed src5"<< norm2(src5)<<" -> "<<norm2(src5_p)<<std::endl;
-
-    Gamma::Algebra Gmu [] = {
-      Gamma::Algebra::GammaX,
-      Gamma::Algebra::GammaY,
-      Gamma::Algebra::GammaZ,
-      Gamma::Algebra::GammaT,
-      Gamma::Algebra::Gamma5
-    };
-    LatticeFermionD    Kinetic(FGrid); Kinetic = Zero();
-    LatticeComplexD    kmu(FGrid);
-    LatticeInteger     scoor(FGrid);
-    LatticeComplexD    sk (FGrid); sk = Zero();
-    LatticeComplexD    sk2(FGrid); sk2= Zero();
-    LatticeComplexD    W(FGrid); W= Zero();
-    LatticeComplexD    one5(FGrid); one5 =ComplexD(1.0,0.0);
-
-    for(int mu=0;mu<Nd;mu++) {
-      LatticeCoordinate(kmu,mu+1);
-      RealD TwoPiL =  M_PI * 2.0/ latt_size[mu];
-      kmu = TwoPiL * kmu;
-      sk2 = sk2 + 2.0*sin(kmu*0.5)*sin(kmu*0.5);
-      sk  = sk  +     sin(kmu)    *sin(kmu);
-      Kinetic = Kinetic + sin(kmu)*ci*(Gamma(Gmu[mu])*src5_p);
-    }
-    std::cout << " src5    "<<norm2(src5_p)<<std::endl;
-    std::cout << " Kinetic "<<norm2(Kinetic)<<std::endl;
-
-    W = one5 - M5 + sk2;
-    std::cout << " W "<<norm2(W)<<std::endl;
-    Kinetic = Kinetic + W * src5_p;
-    std::cout << " Kinetic "<<norm2(Kinetic)<<std::endl;
-
-    LatticeCoordinate(scoor,sdir);
-
-    tmp5 = Cshift(src5_p,sdir,+1);
-    tmp5 = (tmp5 - G5*tmp5)*0.5;
-    tmp5 = where(scoor==Integer(Ls-1),mass*tmp5,-tmp5);
-    Kinetic = Kinetic + tmp5;
-
-    tmp5 = Cshift(src5_p,sdir,-1);
-    tmp5 = (tmp5 + G5*tmp5)*0.5;
-    tmp5 = where(scoor==Integer(0),mass*tmp5,-tmp5);
-    Kinetic = Kinetic + tmp5;
-
-    std::cout<<"Momentum space Ddwf  "<< norm2(Kinetic)<<std::endl;
-    std::cout<<"Stencil Ddwf         "<< norm2(result5)<<std::endl;
-
-    result5 = result5 - Kinetic;
-    std::cout<<"diff "<< norm2(result5)<<std::endl;
-    GRID_ASSERT(norm2(result5)<1.0e-4);
-  }
-
-  ////////////////////////////////////////////////////
-  // Dwf prop
-  ////////////////////////////////////////////////////
-  {
-    std::cout<<"****************************************"<<std::endl;
-    std::cout << "Testing Ddwf Ht Mom space 4d propagator \n";
-    std::cout<<"****************************************"<<std::endl;
-
-    LatticeFermionD    src(&GRID); gaussian(pRNG,src);
-    LatticeFermionD    tmp(&GRID);
-    LatticeFermionD    ref(&GRID);
-    LatticeFermionD    diff(&GRID);
-
-    Coordinate point(4,0);
-    src=Zero();
-    SpinColourVectorD ferm; gaussian(sRNG,ferm);
-    pokeSite(ferm,src,point);
-
-    const int Ls=32;
-    GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,&GRID);
-    GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,&GRID);
-
-    RealD mass=0.01;
-    RealD M5  =0.8;
-    DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,GRID,RBGRID,mass,M5);
-
-    std::cout << " Solving by FFT and Feynman rules" <<std::endl;
-    bool fiveD = false;
-    Ddwf.FreePropagator(src,ref,mass,fiveD);
-
-    Gamma G5(Gamma::Algebra::Gamma5);
-
-    LatticeFermionD    src5(FGrid); src5=Zero();
-    LatticeFermionD    tmp5(FGrid);
-    LatticeFermionD    result5(FGrid); result5=Zero();
-    LatticeFermionD    result4(&GRID);
-    const int sdir=0;
-
-    tmp =   (src + G5*src)*0.5;      InsertSlice(tmp,src5,   0,sdir);
-    tmp =   (src - G5*src)*0.5;      InsertSlice(tmp,src5,Ls-1,sdir);
-
-    std::cout << " Solving by Conjugate Gradient (CGNE)" <<std::endl;
-    Ddwf.Mdag(src5,tmp5);
-    src5=tmp5;
-    MdagMLinearOperator<DomainWallFermionD,LatticeFermionD> HermOp(Ddwf);
-    ConjugateGradient<LatticeFermionD> CG(1.0e-8,10000);
-    CG(HermOp,src5,result5);
-
-    ExtractSlice(tmp,result5,0   ,sdir);   result4 =         (tmp-G5*tmp)*0.5;
-    ExtractSlice(tmp,result5,Ls-1,sdir);   result4 = result4+(tmp+G5*tmp)*0.5;
-
-    std::cout << " Taking difference" <<std::endl;
-    std::cout << "Ddwf result4 "<<norm2(result4)<<std::endl;
-    std::cout << "Ddwf ref     "<<norm2(ref)<<std::endl;
-
-    diff = ref - result4;
-    std::cout << "result - ref     "<<norm2(diff)<<std::endl;
-    GRID_ASSERT(norm2(diff)<1.0e-4);
-  }
-
-  Grid_finalize();
-}

tests/debug/Test_general_coarse.cc

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

tests/debug/Test_reduction.cc

@@ -73,7 +73,36 @@ void testReduction(GridCartesian *grid, GridParallelRNG &rng,
   Field field(grid);
 
   //--------------------------------------------------------------------
-  // a) Timing: Niter timed calls reporting us/call and GB/s.
+  // a) Gaussian random field: sum_gpu (new CUB path) vs sum_gpu_old
+  //    (preserved hand-rolled shared-memory path).  Both promote lanes
+  //    to double internally, so results should agree to near-roundoff.
+  //--------------------------------------------------------------------
+#if defined(GRID_CUDA) || defined(GRID_HIP) || defined(GRID_SYCL)
+  {
+    gaussian(rng, field);
+
+    autoView(v, field, AcceleratorRead);
+    sobj new_result = sum_gpu    (&v[0], osites);
+    sobj old_result = sum_gpu_old(&v[0], osites);
+
+    sobj  diff  = new_result - old_result;
+    RealD diffn = squaredSum(diff);
+    RealD refn  = squaredSum(old_result);
+    RealD reldiff = (refn > 0.0) ? std::sqrt(diffn / refn) : std::sqrt(diffn);
+
+    // Float fields: both paths cast from double to float, expect O(eps_float).
+    // Double fields: ordering differences at most O(V * eps_double).
+    RealD tol = isFloat ? 1e-6 : 1e-10;
+
+    std::cout << GridLogMessage
+              << name << " random reldiff = " << reldiff << std::endl;
+    check(reldiff < tol, name + " random: sum_gpu agrees with sum_gpu_old");
+  }
+#endif
+
+  //--------------------------------------------------------------------
+  // b) Timing: new (CUB/sycl::reduction) vs old (hand-rolled) path.
+  //    Warmup first, then Niter timed calls; report us/call and GB/s.
   //--------------------------------------------------------------------
 #if defined(GRID_CUDA) || defined(GRID_HIP) || defined(GRID_SYCL)
   {
@@ -84,28 +113,38 @@ void testReduction(GridCartesian *grid, GridParallelRNG &rng,
 
     {
       autoView(v, field, AcceleratorRead);
-      for (int i = 0; i < Nwarm; i++) sum_gpu(&v[0], osites);
+      for (int i = 0; i < Nwarm; i++) sum_gpu    (&v[0], osites);
+      for (int i = 0; i < Nwarm; i++) sum_gpu_old(&v[0], osites);
     }
 
-    RealD t_new;
+    RealD t_new, t_old;
     {
       autoView(v, field, AcceleratorRead);
       t_new = -usecond();
       for (int i = 0; i < Niter; i++) sum_gpu(&v[0], osites);
       t_new += usecond();
     }
+    {
+      autoView(v, field, AcceleratorRead);
+      t_old = -usecond();
+      for (int i = 0; i < Niter; i++) sum_gpu_old(&v[0], osites);
+      t_old += usecond();
+    }
 
-    RealD bytes = (RealD)osites * sizeof(vobj);
-    RealD GBs   = bytes / (t_new / Niter) * 1e-3;
+    RealD bytes   = (RealD)osites * sizeof(vobj);
+    RealD GBs_new = bytes / (t_new / Niter) * 1e-3;
+    RealD GBs_old = bytes / (t_old / Niter) * 1e-3;
 
     std::cout << GridLogMessage << name << " timing (" << Niter << " calls):" << std::endl;
     std::cout << GridLogMessage
-              << "  sum_gpu " << t_new/Niter << " us   " << GBs << " GB/s" << std::endl;
+              << "  sum_gpu     " << t_new/Niter << " us   " << GBs_new << " GB/s" << std::endl;
+    std::cout << GridLogMessage
+              << "  sum_gpu_old " << t_old/Niter << " us   " << GBs_old << " GB/s" << std::endl;
   }
 #endif
 
   //--------------------------------------------------------------------
-  // b) Constant field via field = 1.0.
+  // d) Constant field via field = 1.0.
   //
   //    Grid's iMatrix::operator=(scalar) sets only the diagonal, so:
   //      LatticeComplex       -> scalar 1.0        (Ncomp = 1  nonzero per site)

tests/debug/Test_staggered_hdcg.cc

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

tests/debug/Test_staggered_hdcg_ildg.cc

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