Test_extended_meson_field: add view_open timers to measure MemoryManager H2D transfers

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

Grid/algorithms/Algorithms.h

@@ -53,6 +53,7 @@ NAMESPACE_CHECK(approx);
 #include <Grid/algorithms/deflation/MultiRHSBlockCGLinalg.h>
 // Not really deflation, but useful
 #include <Grid/algorithms/blas/MomentumProject.h>
+#include <Grid/algorithms/blas/A2ASpatialSum.h>
 NAMESPACE_CHECK(deflation);
 #include <Grid/algorithms/iterative/ConjugateGradient.h>
 NAMESPACE_CHECK(ConjGrad);

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,17 +65,16 @@ 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);
@@ -83,29 +82,25 @@ 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);
@@ -113,9 +108,7 @@ 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
 
@@ -126,53 +119,45 @@ 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
 
@@ -183,313 +168,325 @@ 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
 
-class FFT {
-private:
-    
-  double flops;
-  double flops_call;
-  uint64_t usec;
-    
-public:
-    
-  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;}    
+struct FFTbase {
+  double        flops;
+  double        flops_call;
+  uint64_t      usec;
+  GridCartesian *_grid;
 
-  FFT ( GridCartesian * grid ) 
-  {
-    flops=0;
-    usec =0;
-  };
-    
-  ~FFT ( void)  {
-    //    delete sgrid;
+  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();
+    }
   }
-    
-  template<class vobj>
-  void FFT_dim_mask(Lattice<vobj> &result,const Lattice<vobj> &source,Coordinate mask,int sign){
+  t_pencil += usecond();
 
-    //    vgrid=result.Grid();
-    //    conformable(result.Grid(),vgrid);
-    //    conformable(source.Grid(),vgrid);
-    const int Ndim = source.Grid()->Nd();
+  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 {
+public:
+  FFT(GridCartesian *grid) : FFTbase(grid) {}
+  ~FFT() {}
+
+  template<class vobj>
+  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;
+    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){
-    const int Ndim = source.Grid()->Nd();
-    Coordinate mask(Ndim,1);
-    FFT_dim_mask(result,source,mask,sign);
+  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);
   }
 
-
   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());
+  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());
 
-    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 vobj::scalar_object sobj;
     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;
 
-    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 */
+    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;
-    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);
 
-    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];
-
-    const Coordinate ldims = grid->_ldimensions;
-    const Coordinate rdims = grid->_rdimensions;
-    const Coordinate sdims = grid->_simd_layout;
-
-    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;
-
-    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);
-
-	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
+    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);
+  }
+};
 
-    t_total +=usecond();
+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::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;
-    
+  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);
+    }
+  }
+
+  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();
+  }
+
+public:
+  PlannedFFT(GridCartesian *grid) : FFTbase(grid) { PlanCreate(); }
+  ~PlannedFFT() { PlanDestroy(); }
+
+  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;
+      }
+    }
+  }
+
+  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_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);
   }
 };
 

Grid/algorithms/blas/A2ASpatialSum.h

@@ -0,0 +1,213 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: Grid/algorithms/blas/A2ASpatialSum.h
+
+    Copyright (C) 2025
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    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
+
+NAMESPACE_BEGIN(Grid);
+
+/*
+  A2ASpatialSum
+
+  Replaces the scalar spatial accumulation loop in A2A extended meson field
+  contractions with a batched GEMM over local time slices, enabling GPU offload.
+
+  Given:
+    leftv[N_i][osite]    — conjugated left SpinColourVectors (SIMD-packed)
+    loopRight[N_j][osite]— type-contracted right SpinColourVectors (SIMD-packed)
+
+  Computes:
+    EMF[i,j,t] = sum_{x,s,c} leftv[i][x,t,s,c] * loopRight[j][x,t,s,c]
+
+  via batched GEMM over nt local time slices, then GlobalSumVector across MPI.
+
+  Memory layout (all C row-major):
+    W_buf  [nt][N_i][nxyz*Nsc]  — W[t][i][x*Nsc+sc]  = leftv[i] at (x,t)
+    LR_buf [nt][N_j][nxyz*Nsc]  — LR[t][j][x*Nsc+sc] = loopRight[j] at (x,t)
+    EMF_buf[nt][N_j][N_i]       — column-major result; EMF[i,j,t] = EMF_buf[t][j][i]
+
+  BLAS call (column-major, OP_T on A so A is read as W[i][k]):
+    C = A^T * B  where A=W[N_i×K C-row], B=LR[N_j×K C-row], C=[N_j×N_i C-row]
+    → C[i,j] = sum_k W[i][k] * LR[j][k] = EMF[i,j] ✓
+*/
+template<class vobj>
+class A2ASpatialSum
+{
+public:
+  typedef typename vobj::scalar_type   scalar;
+  typedef typename vobj::scalar_object sobj;
+
+  GridBase *grid;
+  int N_i, N_j;
+  int nt, nxyz, Nsc;
+
+  deviceVector<scalar>   W_buf;
+  deviceVector<scalar>   LR_buf;
+  deviceVector<scalar>   EMF_buf;
+  deviceVector<scalar *> W_ptrs;
+  deviceVector<scalar *> LR_ptrs;
+  deviceVector<scalar *> EMF_ptrs;
+
+  A2ASpatialSum() : grid(nullptr), N_i(0), N_j(0), nt(0), nxyz(0), Nsc(0) {}
+
+  void Allocate(int _N_i, int _N_j, GridBase *_grid)
+  {
+    grid = _grid;
+    N_i  = _N_i;
+    N_j  = _N_j;
+    Coordinate ldims = grid->LocalDimensions();
+    nt   = ldims[grid->Nd() - 1];
+    nxyz = grid->lSites() / nt;
+    Nsc  = sizeof(sobj) / sizeof(scalar);
+
+    W_buf.resize(nt * N_i * nxyz * Nsc);
+    LR_buf.resize(nt * N_j * nxyz * Nsc);
+    EMF_buf.resize(nt * N_j * N_i);
+
+    // Build persistent batch pointer arrays
+    W_ptrs.resize(nt);
+    LR_ptrs.resize(nt);
+    EMF_ptrs.resize(nt);
+    scalar *Wh   = &W_buf[0];
+    scalar *LRh  = &LR_buf[0];
+    scalar *EMFh = &EMF_buf[0];
+    int lN_i = N_i, lN_j = N_j, lnxyz = nxyz, lNsc = Nsc;
+    for (int t = 0; t < nt; t++) {
+      acceleratorPut(W_ptrs[t],   Wh   + t * lN_i * lnxyz * lNsc);
+      acceleratorPut(LR_ptrs[t],  LRh  + t * lN_j * lnxyz * lNsc);
+      acceleratorPut(EMF_ptrs[t], EMFh + t * lN_j * lN_i);
+    }
+  }
+
+  void PackLeft(const std::vector<Lattice<vobj>> &leftv)
+  {
+    GRID_ASSERT((int)leftv.size() == N_i);
+    PackVectors(leftv, &W_buf[0], N_i);
+  }
+
+  void PackRight(const std::vector<Lattice<vobj>> &loopRight)
+  {
+    GRID_ASSERT((int)loopRight.size() == N_j);
+    PackVectors(loopRight, &LR_buf[0], N_j);
+  }
+
+private:
+  // Pack vecs[N] lattice fields into buf[nt][N][nxyz*Nsc], extracting all SIMD lanes.
+  void PackVectors(const std::vector<Lattice<vobj>> &vecs, scalar *buf, int N)
+  {
+    int nd     = grid->_ndimension;
+    int osites = grid->oSites();
+    int Nsimd  = vobj::Nsimd();
+    int lN     = N;
+    int lNsc   = Nsc;
+    int lnxyz  = nxyz;
+    Coordinate rdimensions = grid->_rdimensions;
+    Coordinate ldims       = grid->LocalDimensions();
+    Coordinate simd        = grid->_simd_layout;
+
+    for (int n = 0; n < N; n++) {
+      autoView(src_v, vecs[n], AcceleratorRead);
+      accelerator_for(sf, osites, Nsimd, {
+#ifdef GRID_SIMT
+        {
+          int lane = acceleratorSIMTlane(Nsimd);
+#else
+          for (int lane = 0; lane < Nsimd; lane++) {
+#endif
+          Coordinate icoor(nd), ocoor(nd), lcoor(nd);
+          Lexicographic::CoorFromIndex(icoor, lane, simd);
+          Lexicographic::CoorFromIndex(ocoor, sf, rdimensions);
+          for (int d = 0; d < nd; d++)
+            lcoor[d] = rdimensions[d] * icoor[d] + ocoor[d];
+
+          int     l_t = lcoor[nd - 1];
+          Coordinate xyz_coor = lcoor;
+          xyz_coor[nd - 1] = 0;
+          int64_t l_xyz;
+          Lexicographic::IndexFromCoor(xyz_coor, l_xyz, ldims);
+
+          sobj    data   = extractLane(lane, src_v[sf]);
+          scalar *data_s = (scalar *)&data;
+
+          int64_t base = (int64_t)l_t * lN * lnxyz * lNsc
+                       + (int64_t)n   * lnxyz * lNsc
+                       + l_xyz * lNsc;
+          for (int sc = 0; sc < lNsc; sc++)
+            buf[base + sc] = data_s[sc];
+        }
+      });
+    }
+  }
+
+public:
+
+  // Batched GEMM + MPI reduction → result[nt_global][N_i][N_j]
+  //
+  // BLAS (column-major, OP_T on A):
+  //   C[N_j×N_i] = A^T[N_i×K] * B[N_j×K]    with K=nxyz*Nsc
+  //   reading A as C row-major [N_i][K] and B as C row-major [N_j][K]
+  //   → C[i,j] = sum_k W[i,k] * LR[j,k] = EMF[i,j] ✓
+  void Sum(Eigen::Tensor<ComplexD, 3> &result)
+  {
+    GridBLAS BLAS;
+
+    int K = nxyz * Nsc;
+    BLAS.gemmBatched(GridBLAS_OP_T, GridBLAS_OP_N,
+                     N_i, N_j, K,
+                     scalar(1.0),
+                     W_ptrs,
+                     LR_ptrs,
+                     scalar(0.0),
+                     EMF_ptrs);
+    BLAS.synchronise();
+
+    // Copy from device and distribute into global-t layout
+    int nt_global = result.dimension(0);
+    int nd        = grid->Nd();
+    int lt_start  = grid->LocalStarts()[nd - 1];
+
+    std::vector<scalar> host_emf(nt * N_j * N_i);
+    acceleratorCopyFromDevice(&EMF_buf[0], host_emf.data(),
+                              nt * N_j * N_i * sizeof(scalar));
+
+    // EMF_buf[t][j*N_i + i] = EMF[i,j] for local t
+    std::vector<scalar> global_emf(nt_global * N_i * N_j, scalar(0.0));
+    for (int lt = 0; lt < nt; lt++) {
+      int gt = lt + lt_start;
+      for (int i = 0; i < N_i; i++)
+      for (int j = 0; j < N_j; j++)
+        global_emf[gt * N_i * N_j + i * N_j + j] = host_emf[lt * N_j * N_i + j * N_i + i];
+    }
+    grid->GlobalSumVector(global_emf.data(), nt_global * N_i * N_j);
+
+    for (int gt = 0; gt < nt_global; gt++)
+    for (int i  = 0; i  < N_i;      i++)
+    for (int j  = 0; j  < N_j;      j++)
+      result(gt, i, j) = global_emf[gt * N_i * N_j + i * N_j + j];
+  }
+};
+
+NAMESPACE_END(Grid);

Grid/lattice/Lattice_reduction.h

@@ -31,9 +31,6 @@ 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,12 +197,15 @@ __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_old(const vobj *lat, Integer osites)
+inline typename vobj::scalar_objectD sumD_gpu_small(const vobj *lat, Integer osites)
 {
   typedef typename vobj::scalar_objectD sobj;
   typedef decltype(lat) Iterator;
-  
+
   Integer nsimd= vobj::Nsimd();
   Integer size = osites*nsimd;
 
@@ -211,46 +214,193 @@ inline typename vobj::scalar_objectD sumD_gpu_small_old(const vobj *lat, Integer
   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;
 }
 
-template <class vobj>
-inline typename vobj::scalar_objectD sumD_gpu_large_old(const vobj *lat, Integer osites)
+// 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)
 {
   typedef typename vobj::vector_type  vector;
   typedef typename vobj::scalar_typeD scalarD;
-  typedef typename vobj::scalar_objectD sobj;
-  sobj ret;
+  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);
   scalarD *ret_p = (scalarD *)&ret;
-  
-  const int words = sizeof(vobj)/sizeof(vector);
 
-  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++) {
+#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
 
-    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_old(const vobj *lat, Integer osites)
+inline typename vobj::scalar_objectD sumD_gpu(const vobj *lat, Integer osites)
 {
   typedef typename vobj::scalar_objectD sobj;
   sobj ret;
@@ -261,9 +411,9 @@ inline typename vobj::scalar_objectD sumD_gpu_old(const vobj *lat, Integer osite
   int ok = getNumBlocksAndThreads(size, sizeof(sobj), numThreads, numBlocks);
   
   if ( ok ) {
-    ret = sumD_gpu_small_old(lat,osites);
+    ret = sumD_gpu_small(lat,osites);
   } else {
-    ret = sumD_gpu_large_old(lat,osites);
+    ret = sumD_gpu_large(lat,osites);
   }
   return ret;
 }
@@ -272,20 +422,20 @@ inline typename vobj::scalar_objectD sumD_gpu_old(const vobj *lat, Integer osite
 // Return as same precision as input performing reduction in double precision though
 /////////////////////////////////////////////////////////////////////////////////////////////////////////
 template <class vobj>
-inline typename vobj::scalar_object sum_gpu_old(const vobj *lat, Integer osites)
+inline typename vobj::scalar_object sum_gpu(const vobj *lat, Integer osites) 
 {
   typedef typename vobj::scalar_object sobj;
   sobj result;
-  result = sumD_gpu_old(lat,osites);
+  result = sumD_gpu(lat,osites);
   return result;
 }
 
 template <class vobj>
-inline typename vobj::scalar_object sum_gpu_large_old(const vobj *lat, Integer osites)
+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_old(lat,osites);
+  result = sumD_gpu_large(lat,osites);
   return result;
 }
 

Grid/lattice/Lattice_reduction_gpu_cub.h

@@ -1,361 +0,0 @@
-/*************************************************************************************
-    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,65 +6,64 @@ NAMESPACE_BEGIN(Grid);
 
 
 template <class vobj>
-inline typename vobj::scalar_objectD sumD_gpu_tensor_old(const vobj *lat, Integer osites)
+inline typename vobj::scalar_objectD sumD_gpu_tensor(const vobj *lat, Integer osites)
 {
-  typedef typename vobj::scalar_object sobj;
+  typedef typename vobj::scalar_object  sobj;
   typedef typename vobj::scalar_objectD sobjD;
 
-  sobj identity; zeroit(identity);
-  sobj ret; zeroit(ret);
-  Integer nsimd= vobj::Nsimd();
-  { 
-    sycl::buffer<sobj, 1> abuff(&ret, {1});
+  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>{osites},
-                      Reduction,
-                      [=] (sycl::id<1> item, auto &sum) {
-                        auto osite   = item[0];
-                        sum +=Reduce(lat[osite]);
-                      });
+      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;
+                       });
     });
   }
-  sobjD dret; convertType(dret,ret);
-  return dret;
+  return ret;
 }
 
 template <class vobj>
-inline typename vobj::scalar_objectD sumD_gpu_large_old(const vobj *lat, Integer osites)
+inline typename vobj::scalar_objectD sumD_gpu_large(const vobj *lat, Integer osites)
 {
-  return sumD_gpu_tensor_old(lat,osites);
+  return sumD_gpu_tensor(lat,osites);
 }
 template <class vobj>
-inline typename vobj::scalar_objectD sumD_gpu_small_old(const vobj *lat, Integer osites)
+inline typename vobj::scalar_objectD sumD_gpu_small(const vobj *lat, Integer osites)
 {
-  return sumD_gpu_large_old(lat,osites);
+  return sumD_gpu_large(lat,osites);
 }
 
 template <class vobj>
-inline typename vobj::scalar_objectD sumD_gpu_old(const vobj *lat, Integer osites)
+inline typename vobj::scalar_objectD sumD_gpu(const vobj *lat, Integer osites)
 {
-  return sumD_gpu_large_old(lat,osites);
+  return sumD_gpu_large(lat,osites);
 }
 
 /////////////////////////////////////////////////////////////////////////////////////////////////////////
 // Return as same precision as input performing reduction in double precision though
 /////////////////////////////////////////////////////////////////////////////////////////////////////////
 template <class vobj>
-inline typename vobj::scalar_object sum_gpu_old(const vobj *lat, Integer osites)
+inline typename vobj::scalar_object sum_gpu(const vobj *lat, Integer osites) 
 {
   typedef typename vobj::scalar_object sobj;
   sobj result;
-  result = sumD_gpu_old(lat,osites);
+  result = sumD_gpu(lat,osites);
   return result;
 }
 
 template <class vobj>
-inline typename vobj::scalar_object sum_gpu_large_old(const vobj *lat, Integer osites)
+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_old(lat,osites);
+  result = sumD_gpu_large(lat,osites);
   return result;
 }
 

Grid/lattice/Lattice_slicesum_core.h

@@ -1,7 +1,6 @@
 #pragma once
 
 #if defined(GRID_CUDA)
-
 #include <cub/cub.cuh>
 #define gpucub cub
 #define gpuError_t cudaError_t
@@ -57,8 +56,13 @@ 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, ::gpucub::Sum(), 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, GRID_CUB_SUM_OP, 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);
@@ -82,11 +86,13 @@ 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,::gpucub::Sum(), zero_init, 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);
   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,6 +113,14 @@ 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/tensors/Tensor_inner.h

@@ -32,6 +32,33 @@ Author: Christoph Lehner <christoph@lhnr.de>
 
 NAMESPACE_BEGIN(Grid);
 
+//////////////////////////////////////
+// innerProductD scalar overloads must be visible before norm2 is defined
+//////////////////////////////////////
+accelerator_inline ComplexD innerProductD(const ComplexF &l,const ComplexF &r){  return innerProduct(l,r); }
+accelerator_inline ComplexD innerProductD(const ComplexD &l,const ComplexD &r){  return innerProduct(l,r); }
+accelerator_inline RealD    innerProductD(const RealD    &l,const RealD    &r){  return innerProduct(l,r); }
+accelerator_inline RealD    innerProductD(const RealF    &l,const RealF    &r){  return innerProduct(l,r); }
+
+accelerator_inline vComplexD innerProductD(const vComplexD &l,const vComplexD &r){  return innerProduct(l,r); }
+accelerator_inline vRealD    innerProductD(const vRealD    &l,const vRealD    &r){  return innerProduct(l,r); }
+accelerator_inline vComplexD innerProductD(const vComplexF &l,const vComplexF &r)
+{
+  vComplexD la,lb;
+  vComplexD ra,rb;
+  Optimization::PrecisionChange::StoD(l.v,la.v,lb.v);
+  Optimization::PrecisionChange::StoD(r.v,ra.v,rb.v);
+  return innerProduct(la,ra) + innerProduct(lb,rb);
+}
+accelerator_inline vRealD innerProductD(const vRealF &l,const vRealF &r)
+{
+  vRealD la,lb;
+  vRealD ra,rb;
+  Optimization::PrecisionChange::StoD(l.v,la.v,lb.v);
+  Optimization::PrecisionChange::StoD(r.v,ra.v,rb.v);
+  return innerProduct(la,ra) + innerProduct(lb,rb);
+}
+
 ///////////////////////////////////////////////////////////////////////////////////////
 // innerProduct Scalar x Scalar -> Scalar
 // innerProduct Vector x Vector -> Scalar
@@ -138,30 +165,8 @@ auto Reduce (const iScalar<l>& lhs) -> iScalar<decltype(Reduce(lhs._internal))>
 
 //////////////////////////////////////
 // innerProductD : if single promote to double and evaluate with sum 2x
+// (scalar/vector overloads are declared above norm2 for ADL visibility)
 //////////////////////////////////////
-accelerator_inline ComplexD innerProductD(const ComplexF &l,const ComplexF &r){  return innerProduct(l,r); }
-accelerator_inline ComplexD innerProductD(const ComplexD &l,const ComplexD &r){  return innerProduct(l,r); }
-accelerator_inline RealD    innerProductD(const RealD    &l,const RealD    &r){  return innerProduct(l,r); }
-accelerator_inline RealD    innerProductD(const RealF    &l,const RealF    &r){  return innerProduct(l,r); }
-
-accelerator_inline vComplexD innerProductD(const vComplexD &l,const vComplexD &r){  return innerProduct(l,r); }
-accelerator_inline vRealD    innerProductD(const vRealD    &l,const vRealD    &r){  return innerProduct(l,r); }
-accelerator_inline vComplexD innerProductD(const vComplexF &l,const vComplexF &r)
-{  
-  vComplexD la,lb;
-  vComplexD ra,rb;
-  Optimization::PrecisionChange::StoD(l.v,la.v,lb.v);
-  Optimization::PrecisionChange::StoD(r.v,ra.v,rb.v);
-  return innerProduct(la,ra) + innerProduct(lb,rb); 
-}
-accelerator_inline vRealD innerProductD(const vRealF &l,const vRealF &r)
-{  
-  vRealD la,lb;
-  vRealD ra,rb;
-  Optimization::PrecisionChange::StoD(l.v,la.v,lb.v);
-  Optimization::PrecisionChange::StoD(r.v,ra.v,rb.v);
-  return innerProduct(la,ra) + innerProduct(lb,rb); 
-}
 
 // Now do it for vector, matrix, scalar
 template<class l,class r,int N> accelerator_inline

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=2;
+uint32_t accelerator_threads=8;
 uint32_t acceleratorThreads(void)       {return accelerator_threads;};
 void     acceleratorThreads(uint32_t t) {accelerator_threads = t;};
 

Grid/threads/Accelerator.h

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

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/Staggered*' `
+STAG_FERMION_FILES=`  find . -name '*.cc' -path '*/instantiation/*' -path '*/instantiation/StaggeredImpl*' `
 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

@@ -0,0 +1,205 @@
+---
+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, 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, mixed-ABI HIP runtime from wrong GTL library (Frontier/ROCm), and deterministic point-to-point reduction trees as a replacement for MPI_Allreduce.
 user-invocable: true
 allowed-tools:
   - Read
@@ -110,6 +110,51 @@ 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:

skills/ref_a2a_emf_work.md

@@ -0,0 +1,61 @@
+---
+name: ref_a2a_emf_work
+description: "A2A Extended Meson Field GPU offload work — status, file locations, pending task"
+metadata: 
+  node_type: memory
+  type: project
+  originSessionId: 956e80aa-401d-481a-80bb-17f8abe1c131
+---
+
+## What was built
+
+`Grid/algorithms/blas/A2ASpatialSum.h` — batched GEMM spatial sum replacing scalar SIMD accumulation. Included via `Grid/algorithms/Algorithms.h`.
+
+`tests/Test_extended_meson_field.cc` — test with class `A2AExtendedMesonFieldRef` containing:
+- CPU reference path (`use_blas=false`)
+- BLAS path (`use_blas=true`) using `A2ASpatialSum`
+- Per-phase timing with `[ref type=N]` / `[blas type=N]` labels
+- 4 contraction types (0-3), all verified at machine precision (~4e-16 rel_err)
+
+## Pending task: GPU offload class
+
+**Goal**: Write `A2AExtendedMesonFieldGPU` in the same test file, replacing all `thread_for` loops with `accelerator_for`-based free function kernels.
+
+The `thread_for` blocks to replace all have the form:
+```cpp
+thread_for(r, rd, {
+  int so = r * grid->_ostride[orthogdim];
+  for (int n = 0; n < e1; n++)
+  for (int b = 0; b < e2; b++) {
+    int ss = so + n * stride + b;
+    // work
+  }
+});
+```
+Replace with `accelerator_for(ss, grid->oSites(), Nsimd, { ... })`.
+
+**Free functions to write** (each takes `Lattice<T>` args, opens views internally):
+- `A2ALoopPropagator` — outerProduct sum (loop build)
+- `A2APackLeftConjugated` — conjugate left fermion fields into `Lattice<SpinColourVector_v>`
+- `A2ALoopLeftContractionType0/1/2/3` — per-site loop × loop propagator → `tloop`
+- `A2ALoopRightContractionType0/1/2/3` — per-site tloop × right → `loopRight[j]`
+
+**Data structure changes required**:
+- `tloopv`: `std::vector<SpinColourMatrix_v>` → `Lattice<SpinColourMatrix_v>` (PropagatorField)
+- `leftv[i]`: `std::vector<SpinColourVector_v>` → `Lattice<SpinColourVector_v>`
+- `loopRight[j]`: `std::vector<SpinColourVector_v>` → `Lattice<SpinColourVector_v>`
+
+**Why**: `std::vector<vobj>` is host memory, not GPU accessible. See [[ref_lattice_vs_vector]].
+
+**`A2ASpatialSum` impact**: `PackLeft`/`PackRight` currently take `std::vector<std::vector<vobj>>`. Once leftv/loopRight become `std::vector<Lattice<vobj>>`, those signatures must change to match.
+
+## Timing on 8.8.8.16 (N_i=N_j=8, Nloop=4, 1 MPI rank)
+
+Dominant costs:
+- `loop_build`: 4-6 ms (outerProduct over 4 propagators)
+- `pack_loopright`: 0.9-2.2 ms (type-dependent)
+- `spatial_sum` (ref): ~1.5 ms
+- `A2ASpatialSum TOTAL`: 2.5-4.3 ms (PackLeft+PackRight dominate GEMM on small volume)
+
+## Related
+[[ref_accelerator_for]] [[ref_coalesced_views]] [[ref_lattice_vs_vector]] [[ref_grid_simt_pattern]]

skills/ref_accelerator_for.md

@@ -0,0 +1,43 @@
+---
+name: ref_accelerator_for
+description: Grid accelerator_for usage — converting block-strided thread_for to GPU-portable oSites loops
+metadata: 
+  node_type: memory
+  type: reference
+  originSessionId: 956e80aa-401d-481a-80bb-17f8abe1c131
+---
+
+## Pattern: block-strided thread_for → accelerator_for over oSites
+
+Old CPU-only pattern (block-strided over orthog dimension):
+```cpp
+thread_for(r, rd, {
+  int so = r * grid->_ostride[orthogdim];
+  for (int n = 0; n < e1; n++)
+  for (int b = 0; b < e2; b++) {
+    int ss = so + n * stride + b;
+    // work on site ss
+  }
+});
+```
+
+GPU-portable replacement:
+```cpp
+accelerator_for(ss, grid->oSites(), Nsimd, {
+  // work on site ss — one SIMT thread per (osite, lane) on GPU
+  // one thread per osite (lane loop implicit via GRID_SIMT) on CPU
+});
+```
+
+Key rules:
+- `accelerator_for(iter, count, Nsimd, body)` — Nsimd is `vobj::Nsimd()` or `grid->Nsimd()`
+- On CPU: expands to `thread_for` over count, `acceleratorSIMTlane` always returns 0 — must use `#ifdef GRID_SIMT` pattern if iterating lanes explicitly (see [[ref_grid_simt_pattern]])
+- On GPU: one SIMT thread per (iter × lane), `acceleratorSIMTlane(Nsimd)` returns actual lane
+- Loop body must capture only scalar/POD by value or via device-accessible pointers; no `std::vector` or host containers inside the body
+- `Coordinate` inside `accelerator_for` must be `AcceleratorVector<int, MaxDims>` (stack-allocated, device-safe) — Grid's `Coordinate` typedef already satisfies this
+
+## Where defined
+`Grid/threads/Accelerator.h` — CPU path ~line 607; GPU paths in conditional blocks above.
+
+## Model file
+`Grid/algorithms/blas/MomentumProject.h` — `ImportVector` is the canonical example of correct `accelerator_for` + SIMD lane extraction.

skills/ref_coalesced_views.md

@@ -0,0 +1,70 @@
+---
+name: ref_coalesced_views
+description: Grid coalescedRead/coalescedWrite and autoView — GPU-portable field access inside accelerator_for
+metadata: 
+  node_type: memory
+  type: reference
+  originSessionId: 956e80aa-401d-481a-80bb-17f8abe1c131
+---
+
+## View access modes
+
+```cpp
+autoView(v, field, AcceleratorRead);   // read-only, device-accessible
+autoView(v, field, AcceleratorWrite);  // write-only, device-accessible
+autoView(v, field, AcceleratorReadWrite); // read-write, device-accessible
+autoView(v, field, CpuRead);           // CPU only (avoids GPU migration)
+autoView(v, field, CpuWrite);          // CPU only
+```
+
+Views must be opened **before** `accelerator_for` and closed (go out of scope) **after**. Never open a view inside the accelerator_for body.
+
+## coalescedRead / coalescedWrite
+
+Inside `accelerator_for(ss, oSites, Nsimd, { ... })`:
+
+```cpp
+auto site = coalescedRead(v[ss]);          // reads SIMT lane; returns scalar_object on GPU, vobj on CPU
+coalescedWrite(v[ss], site);               // writes SIMT lane
+```
+
+- `coalescedRead(v[ss])` calls `v.operator()(ss)` which on GPU returns `extractLane(lane, v[ss])` — one lane per SIMT thread, contiguous across threads → coalesced
+- On CPU returns the full vobj (no lane extraction needed; handled transparently)
+- The returned type is `decltype(coalescedRead(v[ss]))` — use `auto` or match with scalar_object
+
+## Typical kernel pattern
+
+```cpp
+autoView(out_v, out, AcceleratorWrite);
+autoView(in_v,  in,  AcceleratorRead);
+accelerator_for(ss, grid->oSites(), vobj::Nsimd(), {
+  auto x = coalescedRead(in_v[ss]);
+  // modify x ...
+  coalescedWrite(out_v[ss], x);
+});
+```
+
+## Free function kernel signature
+
+```cpp
+template<class vobj>
+void MyKernel(Lattice<vobj> &out, const Lattice<vobj> &in)
+{
+  GridBase *grid = in.Grid();
+  autoView(out_v, out, AcceleratorWrite);
+  autoView(in_v,  in,  AcceleratorRead);
+  accelerator_for(ss, grid->oSites(), vobj::Nsimd(), {
+    auto x = coalescedRead(in_v[ss]);
+    coalescedWrite(out_v[ss], x);
+  });
+}
+```
+
+## What NOT to do
+- Do not access `std::vector` elements inside `accelerator_for` — not device-accessible
+- Do not use `CpuRead`/`CpuWrite` views inside `accelerator_for` — GPU will fault
+- Do not assign to `v[ss]` directly inside `accelerator_for` — use `coalescedWrite`
+- Do not open multiple write views on the same field simultaneously
+
+## Related
+[[ref_accelerator_for]] [[ref_lattice_vs_vector]]

skills/ref_grid_simt_pattern.md

@@ -0,0 +1,47 @@
+---
+name: ref_grid_simt_pattern
+description: Grid GRID_SIMT
+metadata: 
+  node_type: memory
+  type: reference
+  originSessionId: 956e80aa-401d-481a-80bb-17f8abe1c131
+---
+
+## The problem
+
+On CPU, `accelerator_for(sf, oSites, Nsimd, {...})` expands to `thread_for(sf, oSites, {...})` — one thread per osite. `acceleratorSIMTlane(Nsimd)` always returns **0** on CPU. If you need to iterate all Nsimd lanes (e.g. to extract SIMD-packed data), you must loop explicitly on CPU.
+
+On GPU, `accelerator_for` launches one SIMT thread per (osite × lane). `acceleratorSIMTlane(Nsimd)` returns the actual lane index [0, Nsimd).
+
+## Correct pattern (from MomentumProject::ImportVector)
+
+```cpp
+accelerator_for(sf, osites, Nsimd, {
+#ifdef GRID_SIMT
+  {
+    int lane = acceleratorSIMTlane(Nsimd);
+#else
+    for (int lane = 0; lane < Nsimd; lane++) {
+#endif
+      // body using lane
+    }
+  });
+```
+
+- On GPU: `GRID_SIMT` is defined → single-lane body, lane from hardware
+- On CPU: `GRID_SIMT` is not defined → explicit lane loop inside the osite thread
+
+## When is this needed?
+
+Only when you explicitly need the lane index, e.g.:
+- Extracting scalar data from SIMD-packed `vobj` via `extractLane(lane, src[sf])`
+- Computing full local coordinates from (osite, lane) → `Lexicographic::CoorFromIndex(icoor, lane, simd_layout)`
+
+When using `coalescedRead`/`coalescedWrite`, this pattern is **not needed** — those handle lane selection transparently.
+
+## Pitfall that caused a bug
+
+`A2ASpatialSum::PackVectors` originally used `accelerator_for` without the `#ifdef GRID_SIMT` lane loop. On CPU, only lane=0 was extracted, giving wrong norms (~8× too small for `GEN_SIMD_WIDTH=64`, `Nsimd=4`). Fix: add the `#ifdef GRID_SIMT` pattern. See [[ref_accelerator_for]].
+
+## Model file
+`Grid/algorithms/blas/MomentumProject.h`, function `ImportVector`, lines ~166-207.

skills/ref_lattice_vs_vector.md

@@ -0,0 +1,48 @@
+---
+name: ref_lattice_vs_vector
+description: When to use Lattice<T> vs std::vector<T> for GPU-portable field storage in Grid
+metadata: 
+  node_type: memory
+  type: reference
+  originSessionId: 956e80aa-401d-481a-80bb-17f8abe1c131
+---
+
+## Rule
+
+Use `Lattice<vobj>` (or `std::vector<Lattice<vobj>>`) for any field that will be read or written inside `accelerator_for`. `std::vector<vobj>` is host memory and is NOT device-accessible.
+
+## Before vs after GPU offload
+
+```cpp
+// CPU-only (host memory, not GPU accessible)
+std::vector<SpinColourVector_v> tloopv(oSites, Zero());
+// accessed directly: tloopv[ss]
+
+// GPU-portable
+Lattice<SpinColourVector_v> tloop(grid);
+// accessed via view: autoView(tloop_v, tloop, AcceleratorWrite);
+//                    coalescedWrite(tloop_v[ss], val);
+```
+
+## Corollary: function signatures
+
+CPU-only version:
+```cpp
+void PackLeft(const std::vector<std::vector<vobj>> &leftv);
+```
+
+GPU-portable version:
+```cpp
+void PackLeft(const std::vector<Lattice<vobj>> &leftv);
+```
+
+## deviceVector for raw device buffers
+
+`deviceVector<T>` (defined in Grid) is like `std::vector<T>` but in device-accessible memory. Use for raw scalar scratch/pack buffers (e.g. GEMM input/output staging). Not for structured lattice data.
+
+## Pointer arrays for batched BLAS
+
+`deviceVector<scalar *>` holds batch pointer arrays. Populate with `acceleratorPut(ptrs[t], base + offset)` — sets device-side pointer from host. See `A2ASpatialSum::Allocate`.
+
+## Related
+[[ref_coalesced_views]] [[ref_accelerator_for]]

systems/Frontier/config-command

@@ -1,4 +1,3 @@
-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 \
@@ -9,12 +8,13 @@ CLIME=`spack find --paths c-lime@2-3-9 | grep c-lime| cut -c 15-`
 --disable-gparity \
 --disable-fermion-reps \
 --enable-simd=GPU \
---with-gmp=$OLCF_GMP_ROOT \
---with-mpfr=/opt/cray/pe/gcc/mpfr/3.1.4/ \
+--with-gmp=$GMP \
+--with-mpfr=$MPFR \
+--with-openssl=$OPENSSL \
 --disable-fermion-reps \
 CXX=hipcc MPICXX=mpicxx \
-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"
+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"
 
 
 

systems/Frontier/sourceme.sh

@@ -1,28 +1,14 @@
 
 echo spack
-. /autofs/nccs-svm1_home1/paboyle/Crusher/Grid/spack/share/spack/setup-env.sh
+. /autofs/nccs-svm1_home1/paboyle/spack/share/spack/setup-env.sh
 
-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
+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}' `
 
-#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
+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

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=$PREFIX \
-    --with-mpfr=$PREFIX \
+    --with-gmp=$GMP \
+    --with-mpfr=$MPFR \
     --enable-accelerator=cuda \
     --disable-fermion-reps \
     --disable-unified \

systems/Perlmutter/sourceme.sh

@@ -1,4 +1,6 @@
-
 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/11.4
+module load PrgEnv-gnu cpe-cuda cudatoolkit/12.0

systems/mac-arm/config-command-mpi

@@ -3,12 +3,14 @@
 CXX=mpicxx ../../configure \
 	   --enable-simd=GEN \
 	   --enable-comms=mpi-auto \
-	   --enable-Sp=yes \
 	   --enable-unified=yes \
-	   --prefix /Users/peterboyle/QCD/vtk/Grid/install \
+	   --disable-fermion-reps \
+	   --disable-gparity \
+	   --prefix /Users/peterboyle/QCD/Grid-install \
 	   --with-lime=$CLIME \
 	   --with-openssl=$OPENSSL \
 	   --with-gmp=$GMP \
 	   --with-mpfr=$MPFR \
-	   --disable-debug 
+	   --with-fftw=$FFTW \
+	   --disable-debug
 

systems/mac-arm/sourceme.sh

@@ -0,0 +1,7 @@
+source /Users/peterboyle/QCD//Spack/spack//share/spack/setup-env.sh
+export FFTW=`spack find --paths fftw    | grep ^fftw   | 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}' `
+

tests/Test_extended_meson_field.cc

@@ -0,0 +1,841 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: tests/Test_extended_meson_field.cc
+
+Copyright (C) 2015-2025
+
+Author: Peter Boyle <pboyle@bnl.gov>
+Author: Masaaki Tomii <masaaki.tomii@uconn.edu> (original Hadrons kernels)
+
+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
+*************************************************************************************/
+#include "disable_tests_without_instantiations.h"
+#ifdef ENABLE_FERMION_INSTANTIATIONS
+
+#include <Grid/Grid.h>
+#include <Grid/qcd/utils/A2Autils.h>
+
+using namespace Grid;
+
+typedef WilsonImplD FImpl;
+typedef typename FImpl::FermionField FermionField;
+typedef typename FImpl::SiteSpinor   vobj;
+typedef typename vobj::scalar_type   scalar_type;
+typedef typename vobj::vector_type   vector_type;
+typedef iSpinColourMatrix<vector_type> SpinColourMatrix_v;
+typedef iSpinColourVector<vector_type> SpinColourVector_v;
+typedef iSpinMatrix<vector_type>       SpinMatrix_v;
+typedef iSinglet<vector_type>          Scalar_v;
+typedef iSinglet<scalar_type>          Scalar_s;
+typedef Lattice<SpinColourMatrix_v>    PropagatorField;
+
+// CPU reference + optionally batched GEMM spatial sum, ported from
+// Hadrons/Modules/MContraction/A2AExtendedMesonField.hpp
+// (M. Tomii, mtomii/Hadrons:local-2025-edits). Hadrons infrastructure removed.
+// thread_for / CpuRead / orthogdim=3 preserved.
+class A2AExtendedMesonFieldRef
+{
+public:
+  // result is indexed [nt][N_i][N_j].
+  // use_blas=true replaces the scalar spatial accumulation with A2ASpatialSum.
+  static void compute(
+      Eigen::Tensor<ComplexD, 3> &result,
+      const std::vector<FermionField> &left,
+      const std::vector<FermionField> &right,
+      const std::vector<FermionField> &loop1,
+      const std::vector<FermionField> &loop2,
+      const std::vector<Gamma::Algebra> &gamma1,
+      const std::vector<Gamma::Algebra> &gamma2,
+      int type,
+      bool use_blas = false)
+  {
+    GridBase *grid = left[0].Grid();
+
+    const int orthogdim = 3;
+    int rd     = grid->_rdimensions[orthogdim];
+    int ld     = grid->_ldimensions[orthogdim];
+    int Nd     = grid->_ndimension;
+    int Nsimd  = grid->Nsimd();
+
+    int nt  = result.dimension(0);
+    int N_i = (int)left.size();
+    int N_j = (int)right.size();
+
+    std::string tag = std::string(use_blas ? "[blas" : "[ref ") + " type=" + std::to_string(type) + "]";
+    auto Tms = [](double us) { return us * 1e-3; };
+    double t0;
+
+    // ------------------------------------------------------------------
+    // Loop propagator: sum_k outerProduct(loop1[k], loop2[k])
+    // ------------------------------------------------------------------
+    t0 = usecond();
+    PropagatorField loop(grid);
+    loop = Zero();
+    for (unsigned int k = 0; k < loop1.size(); ++k)
+      loop += outerProduct(loop1[k], loop2[k]);
+    std::cout << GridLogMessage << tag << " loop_build:      " << Tms(usecond()-t0) << " ms\n";
+
+    // ------------------------------------------------------------------
+    // Pack conjugated left vectors
+    // ------------------------------------------------------------------
+    t0 = usecond();
+    std::vector<FermionField> leftv(N_i, grid);
+    for (int i = 0; i < N_i; i++)
+      leftv[i] = conjugate(left[i]);
+    std::cout << GridLogMessage << tag << " pack_left:       " << Tms(usecond()-t0) << " ms\n";
+
+    // ------------------------------------------------------------------
+    // Per-site loop contraction into PropagatorField tloop (type-dependent)
+    // ------------------------------------------------------------------
+    t0 = usecond();
+    PropagatorField tloop(grid);
+    tloop = Zero();
+    {
+      autoView(tloopv, tloop, CpuWrite);
+      autoView(loopv,  loop,  CpuRead);
+      if (type == 0) {
+        thread_for(ss, grid->oSites(), {
+          for (int s1 = 0; s1 < Ns; ++s1)
+          for (int s2 = 0; s2 < Ns; ++s2)
+            tloopv[ss]()(s1,s2)(0,0) = loopv[ss]()(s1,s2)(0,0)
+                                     + loopv[ss]()(s1,s2)(1,1)
+                                     + loopv[ss]()(s1,s2)(2,2);
+        });
+      }
+      if (type == 1) {
+        thread_for(ss, grid->oSites(), {
+          tloopv[ss] = Zero();
+          for (int mu = 0; mu < (int)gamma1.size(); ++mu)
+            tloopv[ss] = tloopv[ss] + Gamma(gamma1[mu]) * loopv[ss] * Gamma(gamma2[mu]);
+        });
+      }
+      if (type == 2) {
+        thread_for(ss, grid->oSites(), {
+          tloopv[ss] = Zero();
+          for (int mu = 0; mu < (int)gamma2.size(); ++mu) {
+            SpinColourMatrix_v tmp = Gamma(gamma2[mu]) * loopv[ss];
+            int s1 = mu / Ns;
+            int s2 = mu % Ns;
+            for (int c1 = 0; c1 < Nc; ++c1)
+            for (int c2 = 0; c2 < Nc; ++c2)
+              tloopv[ss]()(s1,s2)(c1,c2) = tmp()(0,0)(c1,c2) + tmp()(1,1)(c1,c2)
+                                          + tmp()(2,2)(c1,c2) + tmp()(3,3)(c1,c2);
+          }
+        });
+      }
+      if (type == 3) {
+        thread_for(ss, grid->oSites(), {
+          SpinMatrix_v spinLoop = Zero();
+          for (int s1 = 0; s1 < Ns; ++s1)
+          for (int s2 = 0; s2 < Ns; ++s2)
+            spinLoop()(s1,s2)() = loopv[ss]()(s1,s2)(0,0)
+                                + loopv[ss]()(s1,s2)(1,1)
+                                + loopv[ss]()(s1,s2)(2,2);
+          tloopv[ss] = Zero();
+          for (int mu = 0; mu < (int)gamma1.size(); ++mu) {
+            SpinMatrix_v tmp2 = Gamma(gamma1[mu]) * spinLoop * Gamma(gamma2[mu]);
+            for (int s1 = 0; s1 < Ns; ++s1)
+            for (int s2 = 0; s2 < Ns; ++s2)
+              tloopv[ss]()(s1,s2)(0,0) = tloopv[ss]()(s1,s2)(0,0) + tmp2()(s1,s2)();
+          }
+        });
+      }
+    }
+    std::cout << GridLogMessage << tag << " tloop:           " << Tms(usecond()-t0) << " ms\n";
+
+    // Select addLoopRight kernel for this type
+    std::function<void(SpinColourVector_v &,
+                       const SpinColourMatrix_v &,
+                       const SpinColourVector_v &,
+                       const std::vector<Gamma::Algebra> &,
+                       const std::vector<Gamma::Algebra> &)> addLoopRight;
+
+    if (type == 0) {
+      addLoopRight = [](SpinColourVector_v &lR,
+                        const SpinColourMatrix_v &loopm,
+                        const SpinColourVector_v &rightv,
+                        const std::vector<Gamma::Algebra> &g1,
+                        const std::vector<Gamma::Algebra> &g2) {
+        SpinMatrix_v spinLoop = Zero();
+        for (int s1 = 0; s1 < Ns; ++s1)
+        for (int s2 = 0; s2 < Ns; ++s2)
+          spinLoop()(s1,s2)() = loopm()(s1,s2)(0,0);
+        for (int mu = 0; mu < (int)g1.size(); ++mu) {
+          SpinMatrix_v GLoop = Gamma(g2[mu]) * spinLoop;
+          auto trGLoop = GLoop()(0,0)() + GLoop()(1,1)() + GLoop()(2,2)() + GLoop()(3,3)();
+          SpinColourVector_v Grightv = Gamma(g1[mu]) * rightv;
+          for (int s = 0; s < Ns; ++s)
+          for (int c = 0; c < Nc; ++c)
+            lR()(s)(c) += Grightv()(s)(c) * trGLoop;
+        }
+      };
+    }
+    if (type == 1) {
+      addLoopRight = [](SpinColourVector_v &lR,
+                        const SpinColourMatrix_v &loopm,
+                        const SpinColourVector_v &rightv,
+                        const std::vector<Gamma::Algebra> &g1,
+                        const std::vector<Gamma::Algebra> &g2) {
+        for (int s = 0; s < Ns; ++s)
+        for (int c = 0; c < Nc; ++c) {
+          lR()(s)(c)
+            += loopm()(s,0)(c,0) * rightv()(0)(0)
+            +  loopm()(s,0)(c,1) * rightv()(0)(1)
+            +  loopm()(s,0)(c,2) * rightv()(0)(2)
+            +  loopm()(s,1)(c,0) * rightv()(1)(0)
+            +  loopm()(s,1)(c,1) * rightv()(1)(1)
+            +  loopm()(s,1)(c,2) * rightv()(1)(2)
+            +  loopm()(s,2)(c,0) * rightv()(2)(0)
+            +  loopm()(s,2)(c,1) * rightv()(2)(1)
+            +  loopm()(s,2)(c,2) * rightv()(2)(2)
+            +  loopm()(s,3)(c,0) * rightv()(3)(0)
+            +  loopm()(s,3)(c,1) * rightv()(3)(1)
+            +  loopm()(s,3)(c,2) * rightv()(3)(2);
+        }
+      };
+    }
+    if (type == 2) {
+      addLoopRight = [](SpinColourVector_v &lR,
+                        const SpinColourMatrix_v &loopm,
+                        const SpinColourVector_v &rightv,
+                        const std::vector<Gamma::Algebra> &g1,
+                        const std::vector<Gamma::Algebra> &g2) {
+        for (int mu = 0; mu < (int)g1.size(); ++mu) {
+          int s1 = mu / Ns;
+          int s2 = mu % Ns;
+          SpinColourVector_v Grightv = Gamma(g1[mu]) * rightv;
+          for (int s = 0; s < Ns; ++s)
+          for (int c = 0; c < Nc; ++c)
+            lR()(s)(c) += loopm()(s1,s2)(c,0) * Grightv()(s)(0)
+                       +  loopm()(s1,s2)(c,1) * Grightv()(s)(1)
+                       +  loopm()(s1,s2)(c,2) * Grightv()(s)(2);
+        }
+      };
+    }
+    if (type == 3) {
+      addLoopRight = [](SpinColourVector_v &lR,
+                        const SpinColourMatrix_v &loopm,
+                        const SpinColourVector_v &rightv,
+                        const std::vector<Gamma::Algebra> &g1,
+                        const std::vector<Gamma::Algebra> &g2) {
+        for (int s = 0; s < Ns; ++s)
+        for (int c = 0; c < Nc; ++c)
+          lR()(s)(c) += loopm()(s,0)(0,0) * rightv()(0)(c)
+                     +  loopm()(s,1)(0,0) * rightv()(1)(c)
+                     +  loopm()(s,2)(0,0) * rightv()(2)(c)
+                     +  loopm()(s,3)(0,0) * rightv()(3)(c);
+      };
+    }
+
+    // ------------------------------------------------------------------
+    // Pack loopRight[j] = type-kernel(tloop, right[j]) per site
+    // ------------------------------------------------------------------
+    t0 = usecond();
+    std::vector<FermionField> loopRight(N_j, grid);
+    {
+      autoView(tlv, tloop, CpuRead);
+      for (int j = 0; j < N_j; j++) {
+        loopRight[j] = Zero();
+        autoView(lRv, loopRight[j], CpuWrite);
+        autoView(rv,  right[j],     CpuRead);
+        thread_for(ss, grid->oSites(), {
+          addLoopRight(lRv[ss], tlv[ss], rv[ss], gamma1, gamma2);
+        });
+      }
+    }
+    std::cout << GridLogMessage << tag << " pack_loopright:  " << Tms(usecond()-t0) << " ms\n";
+
+    if (use_blas) {
+      // ------------------------------------------------------------------
+      // BLAS path: A2ASpatialSum (Allocate + PackLeft + PackRight + Sum)
+      // ------------------------------------------------------------------
+      A2ASpatialSum<SpinColourVector_v> spatial_sum;
+      double t_blas_start = usecond();
+
+      t0 = usecond();
+      spatial_sum.Allocate(N_i, N_j, grid);
+      std::cout << GridLogMessage << tag << " Allocate:        " << Tms(usecond()-t0) << " ms\n";
+
+      t0 = usecond();
+      spatial_sum.PackLeft(leftv);
+      std::cout << GridLogMessage << tag << " PackLeft:        " << Tms(usecond()-t0) << " ms\n";
+
+      t0 = usecond();
+      spatial_sum.PackRight(loopRight);
+      std::cout << GridLogMessage << tag << " PackRight:       " << Tms(usecond()-t0) << " ms\n";
+
+      t0 = usecond();
+      spatial_sum.Sum(result);
+      std::cout << GridLogMessage << tag << " Sum (GEMM+MPI):  " << Tms(usecond()-t0) << " ms\n";
+
+      std::cout << GridLogMessage << tag << " A2ASpatialSum:   " << Tms(usecond()-t_blas_start) << " ms  [TOTAL]\n";
+
+    } else {
+      // ------------------------------------------------------------------
+      // Reference path: SIMD spatial sum + scalar extraction
+      // ------------------------------------------------------------------
+      int MFrvol = rd * N_i * N_j;
+      int MFlvol = ld * N_i * N_j;
+
+      Vector<Scalar_v> lvSum(MFrvol);
+      thread_for(r, MFrvol, { lvSum[r] = Zero(); });
+
+      t0 = usecond();
+      {
+        int e1     = grid->_slice_nblock[orthogdim];
+        int e2     = grid->_slice_block [orthogdim];
+        int stride = grid->_slice_stride[orthogdim];
+        using LView = decltype(leftv[0].View(CpuRead));
+        using RView = decltype(loopRight[0].View(CpuRead));
+        std::vector<LView> lv_views;
+        std::vector<RView> lr_views;
+        lv_views.reserve(N_i);
+        lr_views.reserve(N_j);
+        for (int i = 0; i < N_i; i++) lv_views.push_back(leftv[i].View(CpuRead));
+        for (int j = 0; j < N_j; j++) lr_views.push_back(loopRight[j].View(CpuRead));
+        thread_for(r, rd, {
+          int so   = r * grid->_ostride[orthogdim];
+          int base = N_i * N_j * r;
+          for (int n = 0; n < e1; n++)
+          for (int b = 0; b < e2; b++) {
+            int ss = so + n * stride + b;
+            for (int ii = 0; ii < N_i; ii++)
+            for (int jj = 0; jj < N_j; jj++) {
+              int idx = jj + N_j * ii + base;
+              for (int s = 0; s < Ns; ++s)
+              for (int c = 0; c < Nc; ++c)
+                lvSum[idx]()()() += lv_views[ii][ss]()(s)(c) * lr_views[jj][ss]()(s)(c);
+            }
+          }
+        });
+        for (auto &v : lv_views) v.ViewClose();
+        for (auto &v : lr_views) v.ViewClose();
+      }
+      std::cout << GridLogMessage << tag << " spatial_sum:     " << Tms(usecond()-t0) << " ms\n";
+
+      Vector<Scalar_s> lsSum(MFlvol);
+      thread_for(r, MFlvol, { lsSum[r] = scalar_type(0.0); });
+
+      t0 = usecond();
+      thread_for(rt, rd, {
+        Coordinate icoor(Nd);
+        ExtractBuffer<Scalar_s> extracted(Nsimd);
+        for (int ii = 0; ii < N_i; ii++)
+        for (int jj = 0; jj < N_j; jj++) {
+          int ij_rdx = jj + N_j * (ii + N_i * rt);
+          extract(lvSum[ij_rdx], extracted);
+          for (int idx = 0; idx < Nsimd; idx++) {
+            grid->iCoorFromIindex(icoor, idx);
+            int ldx    = rt + icoor[orthogdim] * rd;
+            int ij_ldx = jj + N_j * (ii + N_i * ldx);
+            lsSum[ij_ldx] = lsSum[ij_ldx] + extracted[idx];
+          }
+        }
+      });
+      std::cout << GridLogMessage << tag << " simd_extract:    " << Tms(usecond()-t0) << " ms\n";
+
+      int pd = grid->_processors[orthogdim];
+      int pc = grid->_processor_coor[orthogdim];
+
+      t0 = usecond();
+      Vector<ComplexD> cache(nt * N_i * N_j, ComplexD(0.0));
+      for (int lt = 0; lt < ld; lt++)
+      for (int pt = 0; pt < pd; pt++) {
+        int t = lt + pt * ld;
+        for (int ii = 0; ii < N_i; ii++)
+        for (int jj = 0; jj < N_j; jj++) {
+          if (pt == pc) {
+            int ij_ldx = jj + N_j * (ii + N_i * lt);
+            cache[jj + N_j * (ii + N_i * t)] = lsSum[ij_ldx]()()();
+          }
+        }
+      }
+      grid->GlobalSumVector(cache.data(), nt * N_i * N_j);
+      std::cout << GridLogMessage << tag << " globalsum:       " << Tms(usecond()-t0) << " ms\n";
+
+      for (int t  = 0; t  < nt;  t++)
+      for (int ii = 0; ii < N_i; ii++)
+      for (int jj = 0; jj < N_j; jj++)
+        result(t, ii, jj) = cache[jj + N_j * (ii + N_i * t)];
+    }
+  }
+};
+
+// ================================================================
+// Free-function GPU kernels — accelerator_for, v(ss) reads,
+// coalescedWrite writes, vobj-level arithmetic throughout.
+// Gamma arrays passed as Vector<Gamma::Algebra> (UVM).
+// ================================================================
+
+void A2ALoopPropagator(PropagatorField &loop,
+                       const std::vector<FermionField> &loop1,
+                       const std::vector<FermionField> &loop2)
+{
+  int Nk      = (int)loop1.size();
+  uint64_t oSites = loop.Grid()->oSites();
+  int Nsimd   = SpinColourVector_v::Nsimd();
+
+  typedef decltype(loop1[0].View(AcceleratorRead)) View;
+  std::vector<View> v1, v2;
+  v1.reserve(Nk); v2.reserve(Nk);
+  for (int k = 0; k < Nk; k++) {
+    v1.push_back(loop1[k].View(AcceleratorRead));
+    v2.push_back(loop2[k].View(AcceleratorRead));
+  }
+
+  deviceVector<SpinColourVector_v *> l1p(Nk), l2p(Nk);
+  for (int k = 0; k < Nk; k++) {
+    acceleratorPut(l1p[k], &v1[k][0]);
+    acceleratorPut(l2p[k], &v2[k][0]);
+  }
+
+  autoView(loopv, loop, AcceleratorWrite);
+  SpinColourVector_v **l1 = &l1p[0];
+  SpinColourVector_v **l2 = &l2p[0];
+  int lNk = Nk;
+  accelerator_for(ss, oSites, Nsimd, {
+    auto result = outerProduct(coalescedRead(l1[0][ss]), coalescedRead(l2[0][ss]));
+    for (int k = 1; k < lNk; k++)
+      result = result + outerProduct(coalescedRead(l1[k][ss]), coalescedRead(l2[k][ss]));
+    coalescedWrite(loopv[ss], result);
+  });
+}
+
+void A2APackLeftConjugated(FermionField &out, const FermionField &in)
+{
+  autoView(outv, out, AcceleratorWrite);
+  autoView(inv,  in,  AcceleratorRead);
+  uint64_t Osites = in.Grid()->oSites();
+  int Nsimd = SpinColourVector_v::Nsimd();
+  accelerator_for(ss, Osites, Nsimd, {
+    coalescedWrite(outv[ss], conjugate(inv(ss)));
+  });
+}
+
+// Type 0: colour-trace stored in (s1,s2)(0,0)
+void A2ALoopLeftContractionType0(PropagatorField &tloop, const PropagatorField &loop)
+{
+  autoView(tloopv, tloop, AcceleratorWrite);
+  autoView(loopv,  loop,  AcceleratorRead);
+  uint64_t Osites = loop.Grid()->oSites();
+  int Nsimd = SpinColourMatrix_v::Nsimd();
+  accelerator_for(ss, Osites, Nsimd, {
+    auto l = loopv(ss);
+    auto tmp = l; tmp = Zero();
+    for (int s1 = 0; s1 < Ns; ++s1)
+    for (int s2 = 0; s2 < Ns; ++s2)
+      tmp()(s1,s2)(0,0) = l()(s1,s2)(0,0) + l()(s1,s2)(1,1) + l()(s1,s2)(2,2);
+    coalescedWrite(tloopv[ss], tmp);
+  });
+}
+
+// Type 1: tloop = sum_mu Gamma(g1[mu]) * loop * Gamma(g2[mu])
+void A2ALoopLeftContractionType1(PropagatorField &tloop, const PropagatorField &loop,
+                                  const Vector<Gamma::Algebra> &gamma1,
+                                  const Vector<Gamma::Algebra> &gamma2)
+{
+  int ng = (int)gamma1.size();
+  const Gamma::Algebra *g1 = gamma1.data();
+  const Gamma::Algebra *g2 = gamma2.data();
+  autoView(tloopv, tloop, AcceleratorWrite);
+  autoView(loopv,  loop,  AcceleratorRead);
+  uint64_t Osites = loop.Grid()->oSites();
+  int Nsimd = SpinColourMatrix_v::Nsimd();
+  accelerator_for(ss, Osites, Nsimd, {
+    auto l = loopv(ss);
+    auto tmp = l; tmp = Zero();
+    for (int mu = 0; mu < ng; ++mu)
+      tmp = tmp + Gamma(g1[mu]) * l * Gamma(g2[mu]);
+    coalescedWrite(tloopv[ss], tmp);
+  });
+}
+
+// Type 2: for mu=[0..ng), s1=mu/Ns, s2=mu%Ns;
+//         tloop(s1,s2)(c1,c2) = Tr_spin( Gamma(g2[mu]) * loop )(c1,c2)
+void A2ALoopLeftContractionType2(PropagatorField &tloop, const PropagatorField &loop,
+                                  const Vector<Gamma::Algebra> &gamma2)
+{
+  int ng = (int)gamma2.size();
+  const Gamma::Algebra *g2 = gamma2.data();
+  autoView(tloopv, tloop, AcceleratorWrite);
+  autoView(loopv,  loop,  AcceleratorRead);
+  uint64_t Osites = loop.Grid()->oSites();
+  int Nsimd = SpinColourMatrix_v::Nsimd();
+  accelerator_for(ss, Osites, Nsimd, {
+    auto l = loopv(ss);
+    auto tmp = l; tmp = Zero();
+    for (int mu = 0; mu < ng; ++mu) {
+      auto gtmp = Gamma(g2[mu]) * l;
+      int s1 = mu / Ns;
+      int s2 = mu % Ns;
+      for (int c1 = 0; c1 < Nc; ++c1)
+      for (int c2 = 0; c2 < Nc; ++c2)
+        tmp()(s1,s2)(c1,c2) = gtmp()(0,0)(c1,c2) + gtmp()(1,1)(c1,c2)
+                            + gtmp()(2,2)(c1,c2) + gtmp()(3,3)(c1,c2);
+    }
+    coalescedWrite(tloopv[ss], tmp);
+  });
+}
+
+// Type 3: colour-trace → spin matrix → sum_mu G1*spinLoop*G2 stored in (s1,s2)(0,0)
+void A2ALoopLeftContractionType3(PropagatorField &tloop, const PropagatorField &loop,
+                                  const Vector<Gamma::Algebra> &gamma1,
+                                  const Vector<Gamma::Algebra> &gamma2)
+{
+  int ng = (int)gamma1.size();
+  const Gamma::Algebra *g1 = gamma1.data();
+  const Gamma::Algebra *g2 = gamma2.data();
+  autoView(tloopv, tloop, AcceleratorWrite);
+  autoView(loopv,  loop,  AcceleratorRead);
+  uint64_t Osites = loop.Grid()->oSites();
+  int Nsimd = SpinColourMatrix_v::Nsimd();
+  accelerator_for(ss, Osites, Nsimd, {
+    typedef decltype(coalescedRead(loopv[0])) calcSCMatrix;
+    typedef iSpinMatrix<typename calcSCMatrix::vector_type> calcSpinMatrix;
+    auto l = loopv(ss);
+    calcSpinMatrix spinLoop; spinLoop = Zero();
+    for (int s1 = 0; s1 < Ns; ++s1)
+    for (int s2 = 0; s2 < Ns; ++s2)
+      spinLoop()(s1,s2)() = l()(s1,s2)(0,0) + l()(s1,s2)(1,1) + l()(s1,s2)(2,2);
+    auto tmp = l; tmp = Zero();
+    for (int mu = 0; mu < ng; ++mu) {
+      calcSpinMatrix tmp2 = Gamma(g1[mu]) * spinLoop * Gamma(g2[mu]);
+      for (int s1 = 0; s1 < Ns; ++s1)
+      for (int s2 = 0; s2 < Ns; ++s2)
+        tmp()(s1,s2)(0,0) = tmp()(s1,s2)(0,0) + tmp2()(s1,s2)();
+    }
+    coalescedWrite(tloopv[ss], tmp);
+  });
+}
+
+// Type 0: loopRight = sum_mu Tr(G2*spinLoop) * G1*right
+//   where spinLoop(s1,s2) = tloop(s1,s2)(0,0)
+void A2ALoopRightContractionType0(FermionField &loopRight,
+                                   const PropagatorField &tloop,
+                                   const FermionField &right,
+                                   const Vector<Gamma::Algebra> &gamma1,
+                                   const Vector<Gamma::Algebra> &gamma2)
+{
+  int ng = (int)gamma1.size();
+  const Gamma::Algebra *g1 = gamma1.data();
+  const Gamma::Algebra *g2 = gamma2.data();
+  autoView(lRv, loopRight, AcceleratorWrite);
+  autoView(tlv, tloop,     AcceleratorRead);
+  autoView(rv,  right,     AcceleratorRead);
+  uint64_t Osites = right.Grid()->oSites();
+  int Nsimd = SpinColourVector_v::Nsimd();
+  accelerator_for(ss, Osites, Nsimd, {
+    typedef decltype(coalescedRead(rv[0]))  calcSCVector;
+    typedef decltype(coalescedRead(tlv[0])) calcSCMatrix;
+    typedef iSpinMatrix<typename calcSCMatrix::vector_type> calcSpinMatrix;
+    auto loopm  = tlv(ss);
+    auto rightv = rv(ss);
+    calcSpinMatrix spinLoop; spinLoop = Zero();
+    for (int s1 = 0; s1 < Ns; ++s1)
+    for (int s2 = 0; s2 < Ns; ++s2)
+      spinLoop()(s1,s2)() = loopm()(s1,s2)(0,0);
+    calcSCVector lR; lR = Zero();
+    for (int mu = 0; mu < ng; ++mu) {
+      auto GLoop   = Gamma(g2[mu]) * spinLoop;
+      auto trGLoop = GLoop()(0,0)() + GLoop()(1,1)() + GLoop()(2,2)() + GLoop()(3,3)();
+      auto Grightv = Gamma(g1[mu]) * rightv;
+      for (int s = 0; s < Ns; ++s)
+      for (int c = 0; c < Nc; ++c)
+        lR()(s)(c) = lR()(s)(c) + Grightv()(s)(c) * trGLoop;
+    }
+    coalescedWrite(lRv[ss], lR);
+  });
+}
+
+// Type 1: loopRight = tloop * right  (SpinColourMatrix * SpinColourVector)
+void A2ALoopRightContractionType1(FermionField &loopRight,
+                                   const PropagatorField &tloop,
+                                   const FermionField &right)
+{
+  autoView(lRv, loopRight, AcceleratorWrite);
+  autoView(tlv, tloop,     AcceleratorRead);
+  autoView(rv,  right,     AcceleratorRead);
+  uint64_t Osites = right.Grid()->oSites();
+  int Nsimd = SpinColourVector_v::Nsimd();
+  accelerator_for(ss, Osites, Nsimd, {
+    coalescedWrite(lRv[ss], tlv(ss) * rv(ss));
+  });
+}
+
+// Type 2: loopRight(s)(c) = sum_{mu,c'} tloop(s1,s2)(c,c') * (G(g1[mu])*right)(s)(c')
+void A2ALoopRightContractionType2(FermionField &loopRight,
+                                   const PropagatorField &tloop,
+                                   const FermionField &right,
+                                   const Vector<Gamma::Algebra> &gamma1)
+{
+  int ng = (int)gamma1.size();
+  const Gamma::Algebra *g1 = gamma1.data();
+  autoView(lRv, loopRight, AcceleratorWrite);
+  autoView(tlv, tloop,     AcceleratorRead);
+  autoView(rv,  right,     AcceleratorRead);
+  uint64_t Osites = right.Grid()->oSites();
+  int Nsimd = SpinColourVector_v::Nsimd();
+  accelerator_for(ss, Osites, Nsimd, {
+    typedef decltype(coalescedRead(rv[0])) calcSCVector;
+    auto loopm  = tlv(ss);
+    auto rightv = rv(ss);
+    calcSCVector lR; lR = Zero();
+    for (int mu = 0; mu < ng; ++mu) {
+      int s1 = mu / Ns;
+      int s2 = mu % Ns;
+      auto Grightv = Gamma(g1[mu]) * rightv;
+      for (int s = 0; s < Ns; ++s)
+      for (int c = 0; c < Nc; ++c)
+        lR()(s)(c) = lR()(s)(c)
+                   + loopm()(s1,s2)(c,0) * Grightv()(s)(0)
+                   + loopm()(s1,s2)(c,1) * Grightv()(s)(1)
+                   + loopm()(s1,s2)(c,2) * Grightv()(s)(2);
+    }
+    coalescedWrite(lRv[ss], lR);
+  });
+}
+
+// Type 3: loopRight(s)(c) = sum_{s'} tloop(s,s')(0,0) * right(s')(c)
+void A2ALoopRightContractionType3(FermionField &loopRight,
+                                   const PropagatorField &tloop,
+                                   const FermionField &right)
+{
+  autoView(lRv, loopRight, AcceleratorWrite);
+  autoView(tlv, tloop,     AcceleratorRead);
+  autoView(rv,  right,     AcceleratorRead);
+  uint64_t Osites = right.Grid()->oSites();
+  int Nsimd = SpinColourVector_v::Nsimd();
+  accelerator_for(ss, Osites, Nsimd, {
+    typedef decltype(coalescedRead(rv[0])) calcSCVector;
+    auto loopm  = tlv(ss);
+    auto rightv = rv(ss);
+    calcSCVector lR; lR = Zero();
+    for (int s = 0; s < Ns; ++s)
+    for (int c = 0; c < Nc; ++c)
+      lR()(s)(c) = loopm()(s,0)(0,0) * rightv()(0)(c)
+                 + loopm()(s,1)(0,0) * rightv()(1)(c)
+                 + loopm()(s,2)(0,0) * rightv()(2)(c)
+                 + loopm()(s,3)(0,0) * rightv()(3)(c);
+    coalescedWrite(lRv[ss], lR);
+  });
+}
+
+// ================================================================
+// GPU-offloaded extended meson field: accelerator_for contractions
+// + A2ASpatialSum GEMM spatial reduction.
+// ================================================================
+class A2AExtendedMesonFieldGPU
+{
+public:
+  static void compute(
+      Eigen::Tensor<ComplexD, 3> &result,
+      const std::vector<FermionField> &left,
+      const std::vector<FermionField> &right,
+      const std::vector<FermionField> &loop1,
+      const std::vector<FermionField> &loop2,
+      const std::vector<Gamma::Algebra> &gamma1_in,
+      const std::vector<Gamma::Algebra> &gamma2_in,
+      int type)
+  {
+    GridBase *grid = left[0].Grid();
+    int N_i = (int)left.size();
+    int N_j = (int)right.size();
+
+    std::string tag = std::string("[gpu  type=") + std::to_string(type) + "]";
+    auto Tms = [](double us) { return us * 1e-3; };
+    double t0;
+
+    Vector<Gamma::Algebra> gamma1(gamma1_in.begin(), gamma1_in.end());
+    Vector<Gamma::Algebra> gamma2(gamma2_in.begin(), gamma2_in.end());
+
+    t0 = usecond();
+    for (auto &f : loop1) { autoView(v, f, AcceleratorRead); }
+    for (auto &f : loop2) { autoView(v, f, AcceleratorRead); }
+    std::cout << GridLogMessage << tag << " view_open_loop:  " << Tms(usecond()-t0) << " ms\n";
+
+    t0 = usecond();
+    PropagatorField loop(grid);
+    A2ALoopPropagator(loop, loop1, loop2);
+    std::cout << GridLogMessage << tag << " loop_build:      " << Tms(usecond()-t0) << " ms\n";
+
+    t0 = usecond();
+    for (int i = 0; i < N_i; i++) { autoView(v, left[i], AcceleratorRead); }
+    std::cout << GridLogMessage << tag << " view_open_left:  " << Tms(usecond()-t0) << " ms\n";
+
+    t0 = usecond();
+    std::vector<FermionField> leftv(N_i, grid);
+    for (int i = 0; i < N_i; i++)
+      A2APackLeftConjugated(leftv[i], left[i]);
+    std::cout << GridLogMessage << tag << " pack_left:       " << Tms(usecond()-t0) << " ms\n";
+
+    t0 = usecond();
+    PropagatorField tloop(grid);
+    tloop = Zero();
+    switch (type) {
+    case 0: A2ALoopLeftContractionType0(tloop, loop);                break;
+    case 1: A2ALoopLeftContractionType1(tloop, loop, gamma1, gamma2); break;
+    case 2: A2ALoopLeftContractionType2(tloop, loop, gamma2);         break;
+    case 3: A2ALoopLeftContractionType3(tloop, loop, gamma1, gamma2); break;
+    }
+    std::cout << GridLogMessage << tag << " tloop:           " << Tms(usecond()-t0) << " ms\n";
+
+    t0 = usecond();
+    { autoView(tlv, tloop, AcceleratorRead); }
+    for (int j = 0; j < N_j; j++) { autoView(rv, right[j], AcceleratorRead); }
+    std::cout << GridLogMessage << tag << " view_open_right: " << Tms(usecond()-t0) << " ms\n";
+
+    t0 = usecond();
+    std::vector<FermionField> loopRight(N_j, grid);
+    for (int j = 0; j < N_j; j++) {
+      switch (type) {
+      case 0: A2ALoopRightContractionType0(loopRight[j], tloop, right[j], gamma1, gamma2); break;
+      case 1: A2ALoopRightContractionType1(loopRight[j], tloop, right[j]);                 break;
+      case 2: A2ALoopRightContractionType2(loopRight[j], tloop, right[j], gamma1);         break;
+      case 3: A2ALoopRightContractionType3(loopRight[j], tloop, right[j]);                 break;
+      }
+    }
+    std::cout << GridLogMessage << tag << " pack_loopright:  " << Tms(usecond()-t0) << " ms\n";
+
+    A2ASpatialSum<SpinColourVector_v> spatial_sum;
+    double t_blas = usecond();
+
+    t0 = usecond();
+    spatial_sum.Allocate(N_i, N_j, grid);
+    std::cout << GridLogMessage << tag << " Allocate:        " << Tms(usecond()-t0) << " ms\n";
+
+    t0 = usecond();
+    spatial_sum.PackLeft(leftv);
+    std::cout << GridLogMessage << tag << " PackLeft:        " << Tms(usecond()-t0) << " ms\n";
+
+    t0 = usecond();
+    spatial_sum.PackRight(loopRight);
+    std::cout << GridLogMessage << tag << " PackRight:       " << Tms(usecond()-t0) << " ms\n";
+
+    t0 = usecond();
+    spatial_sum.Sum(result);
+    std::cout << GridLogMessage << tag << " Sum (GEMM+MPI):  " << Tms(usecond()-t0) << " ms\n";
+
+    std::cout << GridLogMessage << tag << " A2ASpatialSum:   " << Tms(usecond()-t_blas) << " ms  [TOTAL]\n";
+  }
+};
+
+int main(int argc, char *argv[])
+{
+  Grid_init(&argc, &argv);
+
+  Coordinate latt_size   = GridDefaultLatt();
+  Coordinate simd_layout = GridDefaultSimd(4, vComplexD::Nsimd());
+  Coordinate mpi_layout  = GridDefaultMpi();
+  GridCartesian grid(latt_size, simd_layout, mpi_layout);
+
+  int Nt    = latt_size[Tp];
+  int N_i   = 8;
+  int N_j   = 8;
+  int Nloop = 4;
+
+  if (GridCmdOptionExists(argv, argv+argc, "--Ni"))
+    N_i = std::stoi(GridCmdOptionPayload(argv, argv+argc, "--Ni"));
+  if (GridCmdOptionExists(argv, argv+argc, "--Nj"))
+    N_j = std::stoi(GridCmdOptionPayload(argv, argv+argc, "--Nj"));
+
+  GridParallelRNG pRNG(&grid);
+  pRNG.SeedFixedIntegers({1, 2, 3, 4});
+
+  std::vector<FermionField> left(N_i,   &grid);
+  std::vector<FermionField> right(N_j,  &grid);
+  std::vector<FermionField> loop1(Nloop, &grid);
+  std::vector<FermionField> loop2(Nloop, &grid);
+
+  for (auto &f : left)  random(pRNG, f);
+  for (auto &f : right) random(pRNG, f);
+  for (auto &f : loop1) random(pRNG, f);
+  for (auto &f : loop2) random(pRNG, f);
+
+  std::vector<Gamma::Algebra> GammaMU = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
+  };
+
+  Eigen::Tensor<ComplexD, 3> result_ref(Nt, N_i, N_j);
+  Eigen::Tensor<ComplexD, 3> result_blas(Nt, N_i, N_j);
+  Eigen::Tensor<ComplexD, 3> result_gpu(Nt, N_i, N_j);
+  double t_ref = 0, t_blas = 0, t_gpu = 0, start, stop;
+
+  for (int type = 0; type < 4; type++) {
+
+    result_ref.setZero();
+    start = usecond();
+    A2AExtendedMesonFieldRef::compute(result_ref,  left, right, loop1, loop2,
+                                      GammaMU, GammaMU, type, false);
+    stop = usecond(); t_ref = stop - start;
+
+    result_blas.setZero();
+    start = usecond();
+    A2AExtendedMesonFieldRef::compute(result_blas, left, right, loop1, loop2,
+                                      GammaMU, GammaMU, type, true);
+    stop = usecond(); t_blas = stop - start;
+
+    result_gpu.setZero();
+    start = usecond();
+    A2AExtendedMesonFieldGPU::compute(result_gpu, left, right, loop1, loop2,
+                                      GammaMU, GammaMU, type);
+    stop = usecond(); t_gpu = stop - start;
+
+    double norm2_ref = 0.0, norm2_blas = 0.0, norm2_gpu = 0.0;
+    double norm2_diff_blas = 0.0, norm2_diff_gpu = 0.0;
+    for (int t  = 0; t  < Nt;  t++)
+    for (int ii = 0; ii < N_i; ii++)
+    for (int jj = 0; jj < N_j; jj++) {
+      norm2_ref  += norm2(result_ref(t, ii, jj));
+      norm2_blas += norm2(result_blas(t, ii, jj));
+      norm2_gpu  += norm2(result_gpu(t, ii, jj));
+      ComplexD diff_blas = result_ref(t, ii, jj) - result_blas(t, ii, jj);
+      ComplexD diff_gpu  = result_ref(t, ii, jj) - result_gpu(t, ii, jj);
+      norm2_diff_blas += norm2(diff_blas);
+      norm2_diff_gpu  += norm2(diff_gpu);
+    }
+
+    double rel_blas = (norm2_ref > 0) ? std::sqrt(norm2_diff_blas / norm2_ref) : 0.0;
+    double rel_gpu  = (norm2_ref > 0) ? std::sqrt(norm2_diff_gpu  / norm2_ref) : 0.0;
+
+    std::cout << GridLogMessage
+              << "type=" << type
+              << "  norm2_ref="  << norm2_ref
+              << "  norm2_blas=" << norm2_blas
+              << "  norm2_gpu="  << norm2_gpu
+              << "  rel_blas="   << rel_blas
+              << "  rel_gpu="    << rel_gpu
+              << "  t_ref="  << t_ref  * 1e-6 << "s"
+              << "  t_blas=" << t_blas * 1e-6 << "s"
+              << "  t_gpu="  << t_gpu  * 1e-6 << "s"
+              << std::endl;
+
+    GRID_ASSERT(rel_blas < 1e-10);
+    GRID_ASSERT(rel_gpu  < 1e-10);
+  }
+
+  std::cout << GridLogMessage << "All types passed A2ASpatialSum and GPU regression." << std::endl;
+
+  Grid_finalize();
+  return EXIT_SUCCESS;
+}
+
+#endif

tests/core/Test_planned_fft.cc

@@ -0,0 +1,321 @@
+/*************************************************************************************
+
+    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_hipfft_bug_fail.cc

@@ -0,0 +1,76 @@
+/*
+ * Isolating the hipfft HIPFFT_PARSE_ERROR on ROCm 7 / hipFFT 1.0.20.
+ *
+ * Tests three orderings with an empty rocFFT cache to find which GPU
+ * operation before plan creation triggers the failure:
+ *   A) hipMalloc only            — hypothesis: passes (no async GPU work)
+ *   B) hipMalloc + hipMemset     — hypothesis: fails  (async GPU work in flight)
+ *   C) hipMalloc + hipMemset     — hypothesis: passes (work completed before plan)
+ *      + hipDeviceSynchronize
+ *
+ * Compile:
+ *   hipcc -o Test_hipfft_bug_fail Test_hipfft_bug_fail.cc -lhipfft
+ *
+ * Run with empty cache:
+ *   rm -rf ~/.cache/
+ *   ./Test_hipfft_bug_fail
+ */
+
+#include <cstdio>
+#include <hipfft/hipfft.h>
+#include <hip/hip_runtime.h>
+
+static const char *res(hipfftResult rv) {
+  return rv == HIPFFT_SUCCESS ? "SUCCESS" : "PARSE_ERROR";
+}
+
+static hipfftResult makePlan(int G, int howmany) {
+  int n[] = {G};
+  hipfftHandle p;
+  size_t workSize = 0;
+  hipfftCreate(&p);
+  hipfftResult rv = hipfftMakePlanMany(p, 1, n,
+                                       nullptr, 1, G, nullptr, 1, G,
+                                       HIPFFT_Z2Z, howmany, &workSize);
+  hipfftDestroy(p);
+  return rv;
+}
+
+int main(void) {
+  hipDeviceProp_t prop;
+  hipGetDeviceProperties(&prop, 0);
+  printf("Device: %s\n", prop.name);
+#ifdef hipfftVersionMinor
+  printf("hipFFT version: %d.%d.%d\n\n",
+         hipfftVersionMajor, hipfftVersionMinor, hipfftVersionPatch);
+#endif
+
+  for (int G : {4, 8, 16, 32}) {
+    int howmany = 512;
+    long nelems = (long)G * howmany;
+    hipfftDoubleComplex *buf = nullptr;
+    hipMalloc(&buf, nelems * sizeof(hipfftDoubleComplex));
+
+    // Tests ordered so each runs before a prior success can populate the cache.
+
+    // B first: hipMalloc + hipMemset (async GPU work in flight)
+    // If this fails, A (no hipMemset) will pass, confirming hipMemset is the trigger.
+    hipMemset(buf, 0, nelems * sizeof(hipfftDoubleComplex));
+    hipfftResult rvB = makePlan(G, howmany);
+    printf("G=%-4d  B) hipMalloc + hipMemset       : %s\n", G, res(rvB));
+
+    // C: hipMalloc + hipMemset + sync — does syncing before plan creation fix it?
+    hipMemset(buf, 0, nelems * sizeof(hipfftDoubleComplex));
+    hipDeviceSynchronize();
+    hipfftResult rvC = makePlan(G, howmany);
+    printf("G=%-4d  C) hipMalloc + hipMemset + sync: %s\n", G, res(rvC));
+
+    // A last: hipMalloc only, no async GPU work — should always pass
+    hipfftResult rvA = makePlan(G, howmany);
+    printf("G=%-4d  A) hipMalloc only             : %s\n\n", G, res(rvA));
+
+    hipFree(buf);
+  }
+
+  return 0;
+}

tests/debug/Test_hipfft_bug_pass.cc

@@ -0,0 +1,61 @@
+/*
+ * Minimal program demonstrating the workaround for the hipfft ROCm 7 bug.
+ *
+ * Workaround: create the hipfft plan BEFORE any hipMalloc.  Plan creation
+ * for G < 32 then succeeds even with an empty rocFFT cache.
+ *
+ * Compile:
+ *   hipcc -o Test_hipfft_bug_pass Test_hipfft_bug_pass.cc -lhipfft
+ *
+ * Run:
+ *   rm -rf ~/.cache/rocfft
+ *   ./Test_hipfft_bug_pass
+ *
+ * Expected: all G values succeed.
+ * Compare with Test_hipfft_bug_fail.cc which uses the opposite ordering.
+ */
+
+#include <cstdio>
+#include <hipfft/hipfft.h>
+#include <hip/hip_runtime.h>
+
+int main(void) {
+  hipDeviceProp_t prop;
+  hipGetDeviceProperties(&prop, 0);
+  printf("Device: %s\n", prop.name);
+#ifdef hipfftVersionMinor
+  printf("hipFFT version: %d.%d.%d\n\n",
+         hipfftVersionMajor, hipfftVersionMinor, hipfftVersionPatch);
+#endif
+
+  for (int G : {8, 16, 32}) {
+    int howmany = 512;
+    int n[] = {G};
+    long nelems = (long)G * howmany;
+
+    // Plan created BEFORE hipMalloc — succeeds for all G
+    hipfftHandle p;
+    size_t workSize = 0;
+    hipfftCreate(&p);
+    hipfftResult rv = hipfftMakePlanMany(p, 1, n,
+                                         nullptr, 1, G, nullptr, 1, G,
+                                         HIPFFT_Z2Z, howmany, &workSize);
+    printf("G=%-4d  plan-then-hipMalloc: %d (%s)\n",
+           G, (int)rv, rv == HIPFFT_SUCCESS ? "HIPFFT_SUCCESS" : "HIPFFT_PARSE_ERROR");
+
+    if (rv == HIPFFT_SUCCESS) {
+      hipfftDoubleComplex *buf = nullptr;
+      hipMalloc(&buf, nelems * sizeof(hipfftDoubleComplex));
+      hipMemset(buf, 0, nelems * sizeof(hipfftDoubleComplex));
+      rv = hipfftExecZ2Z(p, buf, buf, HIPFFT_FORWARD);
+      hipDeviceSynchronize();
+      printf("G=%-4d  execFwd:             %d (%s)\n",
+             G, (int)rv, rv == HIPFFT_SUCCESS ? "HIPFFT_SUCCESS" : "FAILED");
+      hipFree(buf);
+    }
+
+    hipfftDestroy(p);
+  }
+
+  return 0;
+}

tests/debug/Test_hipfft_minimal.cc

@@ -0,0 +1,161 @@
+/*
+ * Minimal reproducer for hipfftMakePlanMany / hipfftPlanMany failures.
+ *
+ * Compile on Frontier (no Grid headers needed):
+ *   hipcc -o Test_hipfft_minimal Test_hipfft_minimal.cc -lhipfft
+ *
+ * Run:
+ *   ./Test_hipfft_minimal
+ */
+
+#include <cstdio>
+#include <cstdlib>
+#include <hipfft/hipfft.h>
+#include <hip/hip_runtime.h>
+
+static const char *hipfftResultString(hipfftResult r) {
+  switch (r) {
+  case HIPFFT_SUCCESS:                  return "HIPFFT_SUCCESS";
+  case HIPFFT_INVALID_PLAN:             return "HIPFFT_INVALID_PLAN";
+  case HIPFFT_ALLOC_FAILED:             return "HIPFFT_ALLOC_FAILED";
+  case HIPFFT_INVALID_TYPE:             return "HIPFFT_INVALID_TYPE";
+  case HIPFFT_INVALID_VALUE:            return "HIPFFT_INVALID_VALUE";
+  case HIPFFT_INTERNAL_ERROR:           return "HIPFFT_INTERNAL_ERROR";
+  case HIPFFT_EXEC_FAILED:              return "HIPFFT_EXEC_FAILED";
+  case HIPFFT_SETUP_FAILED:             return "HIPFFT_SETUP_FAILED";
+  case HIPFFT_INVALID_SIZE:             return "HIPFFT_INVALID_SIZE";
+  case HIPFFT_UNALIGNED_DATA:           return "HIPFFT_UNALIGNED_DATA";
+  case HIPFFT_INCOMPLETE_PARAMETER_LIST:return "HIPFFT_INCOMPLETE_PARAMETER_LIST";
+  case HIPFFT_INVALID_DEVICE:           return "HIPFFT_INVALID_DEVICE";
+  case HIPFFT_PARSE_ERROR:              return "HIPFFT_PARSE_ERROR";
+  case HIPFFT_NO_WORKSPACE:             return "HIPFFT_NO_WORKSPACE";
+  case HIPFFT_NOT_IMPLEMENTED:          return "HIPFFT_NOT_IMPLEMENTED";
+  case HIPFFT_NOT_SUPPORTED:            return "HIPFFT_NOT_SUPPORTED";
+  default:                              return "UNKNOWN";
+  }
+}
+
+// Plan creation + execution for (G, howmany).
+// Tests two orderings to isolate whether a prior hipMalloc poisons hipfft
+// plan creation for small G on ROCm 7:
+//   A) plan BEFORE hipMalloc  — hypothesis: succeeds
+//   B) hipMalloc BEFORE plan  — hypothesis: fails for G < 32
+static void tryPlanAndExec(int G, long howmany) {
+  int n[] = {G};
+  long nelems = (long)G * howmany;
+
+  printf("--- G=%-4d  howmany=%-10ld  total_elems=%-12ld ---\n",
+         G, howmany, nelems);
+
+  // Allocate device buffer (hipfftDoubleComplex = 16 bytes each)
+  hipfftDoubleComplex *dbuf = nullptr;
+  hipError_t herr = hipMalloc(&dbuf, nelems * sizeof(hipfftDoubleComplex));
+  if (herr != hipSuccess) {
+    printf("  hipMalloc failed (%d) for %ld elems — skipping\n\n", (int)herr, nelems);
+    return;
+  }
+  hipMemset(dbuf, 0, nelems * sizeof(hipfftDoubleComplex));
+
+  // 1. hipfftPlanMany (one-step, nullptr embed) — current Grid path
+  {
+    hipfftHandle p;
+    hipfftResult rv = hipfftPlanMany(&p, 1, n,
+                                     nullptr, 1, G,
+                                     nullptr, 1, G,
+                                     HIPFFT_Z2Z, (int)howmany);
+    printf("  hipfftPlanMany   create : %d (%s)\n", (int)rv, hipfftResultString(rv));
+    if (rv == HIPFFT_SUCCESS) {
+      rv = hipfftExecZ2Z(p, dbuf, dbuf, HIPFFT_FORWARD);
+      hipDeviceSynchronize();
+      printf("  hipfftPlanMany   execFwd: %d (%s)\n", (int)rv, hipfftResultString(rv));
+      hipfftDestroy(p);
+    }
+  }
+
+  // 2. hipfftCreate + hipfftMakePlanMany (two-step) — also current Grid path
+  {
+    hipfftHandle p;
+    size_t workSize = 0;
+    hipfftResult rc = hipfftCreate(&p);
+    if (rc == HIPFFT_SUCCESS) {
+      hipfftResult rv = hipfftMakePlanMany(p, 1, n,
+                                           nullptr, 1, G,
+                                           nullptr, 1, G,
+                                           HIPFFT_Z2Z, (int)howmany, &workSize);
+      printf("  hipfftMakePlanMany      : %d (%s)  workSize=%zu\n",
+             (int)rv, hipfftResultString(rv), workSize);
+      if (rv == HIPFFT_SUCCESS) {
+        rv = hipfftExecZ2Z(p, dbuf, dbuf, HIPFFT_FORWARD);
+        hipDeviceSynchronize();
+        printf("  hipfftMakePlanMany exec : %d (%s)\n", (int)rv, hipfftResultString(rv));
+      }
+      hipfftDestroy(p);
+    } else {
+      printf("  hipfftCreate            : %d (%s)\n", (int)rc, hipfftResultString(rc));
+    }
+  }
+
+  // 3. hipfftPlan1d (simplest API, batch = howmany)
+  {
+    hipfftHandle p;
+    hipfftResult rv = hipfftPlan1d(&p, G, HIPFFT_Z2Z, (int)howmany);
+    printf("  hipfftPlan1d     create : %d (%s)\n", (int)rv, hipfftResultString(rv));
+    if (rv == HIPFFT_SUCCESS) {
+      rv = hipfftExecZ2Z(p, dbuf, dbuf, HIPFFT_FORWARD);
+      hipDeviceSynchronize();
+      printf("  hipfftPlan1d     execFwd: %d (%s)\n", (int)rv, hipfftResultString(rv));
+      hipfftDestroy(p);
+    }
+  }
+
+  hipFree(dbuf);
+  printf("\n");
+}
+
+int main(void) {
+  // Print HIP device info
+  int device = 0;
+  hipGetDevice(&device);
+  hipDeviceProp_t prop;
+  hipGetDeviceProperties(&prop, device);
+  printf("Device %d: %s  warpSize=%d\n\n", device, prop.name, prop.warpSize);
+
+#ifdef hipfftVersionMinor
+  printf("hipFFT version: %d.%d.%d\n\n",
+         hipfftVersionMajor, hipfftVersionMinor, hipfftVersionPatch);
+#endif
+
+  // Original sweep with small howmany (these passed first time)
+  printf("=== Small howmany (original sweep) ===\n\n");
+  for (int G : {4, 8, 12, 16, 24, 32, 48, 64})
+    tryPlanAndExec(G, 512);
+
+  // Grid-realistic howmany values derived from actual lattice geometries.
+  // howmany = Ncomp * product(ldimensions[d] for d != dim)
+  // For LatticeComplexD: Ncomp=1.
+  printf("=== Grid-realistic parameters ===\n\n");
+
+  // --grid 16.16.16.16  4D FFT  (KNOWN TO FAIL in Grid)
+  // Each dim: G=16, Nperp=16^3=4096
+  tryPlanAndExec(16, 4096);
+
+  // --grid 32.32.32.32  4D FFT  (KNOWN TO SUCCEED in Grid)
+  // Each dim: G=32, Nperp=32^3=32768
+  tryPlanAndExec(32, 32768);
+
+  // --grid 32.32.32.32 Ls=8  5D DWF FFT  (KNOWN TO FAIL on dim 0 in Grid)
+  // dim 0: G=8,  Nperp=32^4=1048576
+  tryPlanAndExec(8, 1048576);
+  // dim 1-4: G=32, Nperp=8*32^3=262144
+  tryPlanAndExec(32, 262144);
+
+  // Extra intermediate cases to bracket the failure
+  tryPlanAndExec(16, 1024);
+  tryPlanAndExec(16, 2048);
+  tryPlanAndExec(16, 8192);
+  tryPlanAndExec(8,  4096);
+  tryPlanAndExec(8,  65536);
+  tryPlanAndExec(8,  262144);
+
+  return 0;
+}

tests/debug/Test_hipfft_repro.cc

@@ -0,0 +1,168 @@
+/*
+ * Reproducer for HIPFFT_PARSE_ERROR (error 12) from hipfftMakePlanMany on
+ * ROCm 7 / hipFFT 1.0.20 (Frontier, MI210 login and MI250X compute nodes).
+ *
+ * Observed failure: G < 32 returns HIPFFT_PARSE_ERROR from all three plan
+ * creation APIs (hipfftPlanMany, hipfftMakePlanMany, hipfftPlan1d) when a
+ * device buffer is allocated and zeroed with hipMalloc+hipMemset before the
+ * plan creation call.  G >= 32 succeeds.
+ *
+ * Contrast with Test_hipfft_minimal.cc (plan-first ordering) which passes
+ * for all G even with an empty rocFFT cache.
+ *
+ * Compile on Frontier (no Grid headers needed):
+ *   hipcc -o Test_hipfft_repro Test_hipfft_repro.cc -lhipfft
+ *
+ * Run with empty cache to reproduce the failure:
+ *   rm -rf ~/.cache/rocfft
+ *   ./Test_hipfft_repro
+ */
+
+#include <cstdio>
+#include <cstdlib>
+#include <hipfft/hipfft.h>
+#include <hip/hip_runtime.h>
+
+static const char *hipfftResultString(hipfftResult r) {
+  switch (r) {
+  case HIPFFT_SUCCESS:                  return "HIPFFT_SUCCESS";
+  case HIPFFT_INVALID_PLAN:             return "HIPFFT_INVALID_PLAN";
+  case HIPFFT_ALLOC_FAILED:             return "HIPFFT_ALLOC_FAILED";
+  case HIPFFT_INVALID_TYPE:             return "HIPFFT_INVALID_TYPE";
+  case HIPFFT_INVALID_VALUE:            return "HIPFFT_INVALID_VALUE";
+  case HIPFFT_INTERNAL_ERROR:           return "HIPFFT_INTERNAL_ERROR";
+  case HIPFFT_EXEC_FAILED:              return "HIPFFT_EXEC_FAILED";
+  case HIPFFT_SETUP_FAILED:             return "HIPFFT_SETUP_FAILED";
+  case HIPFFT_INVALID_SIZE:             return "HIPFFT_INVALID_SIZE";
+  case HIPFFT_UNALIGNED_DATA:           return "HIPFFT_UNALIGNED_DATA";
+  case HIPFFT_INCOMPLETE_PARAMETER_LIST:return "HIPFFT_INCOMPLETE_PARAMETER_LIST";
+  case HIPFFT_INVALID_DEVICE:           return "HIPFFT_INVALID_DEVICE";
+  case HIPFFT_PARSE_ERROR:              return "HIPFFT_PARSE_ERROR";
+  case HIPFFT_NO_WORKSPACE:             return "HIPFFT_NO_WORKSPACE";
+  case HIPFFT_NOT_IMPLEMENTED:          return "HIPFFT_NOT_IMPLEMENTED";
+  case HIPFFT_NOT_SUPPORTED:            return "HIPFFT_NOT_SUPPORTED";
+  default:                              return "UNKNOWN";
+  }
+}
+
+// Plan creation + execution for (G, howmany) using hipfftCreate+hipfftMakePlanMany.
+// This is the path Grid's FFT.h now uses.
+static void tryPlanAndExec(int G, long howmany) {
+  int n[] = {G};
+  long nelems = (long)G * howmany;
+
+  printf("--- G=%-4d  howmany=%-10ld  total_elems=%-12ld ---\n",
+         G, howmany, nelems);
+
+  // Allocate device buffer (hipfftDoubleComplex = 16 bytes each)
+  hipfftDoubleComplex *dbuf = nullptr;
+  hipError_t herr = hipMalloc(&dbuf, nelems * sizeof(hipfftDoubleComplex));
+  if (herr != hipSuccess) {
+    printf("  hipMalloc failed (%d) for %ld elems — skipping\n\n", (int)herr, nelems);
+    return;
+  }
+  hipMemset(dbuf, 0, nelems * sizeof(hipfftDoubleComplex));
+
+  // 1. hipfftPlanMany (one-step, nullptr embed) — current Grid path
+  {
+    hipfftHandle p;
+    hipfftResult rv = hipfftPlanMany(&p, 1, n,
+                                     nullptr, 1, G,
+                                     nullptr, 1, G,
+                                     HIPFFT_Z2Z, (int)howmany);
+    printf("  hipfftPlanMany   create : %d (%s)\n", (int)rv, hipfftResultString(rv));
+    if (rv == HIPFFT_SUCCESS) {
+      rv = hipfftExecZ2Z(p, dbuf, dbuf, HIPFFT_FORWARD);
+      hipDeviceSynchronize();
+      printf("  hipfftPlanMany   execFwd: %d (%s)\n", (int)rv, hipfftResultString(rv));
+      hipfftDestroy(p);
+    }
+  }
+
+  // 2. hipfftCreate + hipfftMakePlanMany (two-step) — also current Grid path
+  {
+    hipfftHandle p;
+    size_t workSize = 0;
+    hipfftResult rc = hipfftCreate(&p);
+    if (rc == HIPFFT_SUCCESS) {
+      hipfftResult rv = hipfftMakePlanMany(p, 1, n,
+                                           nullptr, 1, G,
+                                           nullptr, 1, G,
+                                           HIPFFT_Z2Z, (int)howmany, &workSize);
+      printf("  hipfftMakePlanMany      : %d (%s)  workSize=%zu\n",
+             (int)rv, hipfftResultString(rv), workSize);
+      if (rv == HIPFFT_SUCCESS) {
+        rv = hipfftExecZ2Z(p, dbuf, dbuf, HIPFFT_FORWARD);
+        hipDeviceSynchronize();
+        printf("  hipfftMakePlanMany exec : %d (%s)\n", (int)rv, hipfftResultString(rv));
+      }
+      hipfftDestroy(p);
+    } else {
+      printf("  hipfftCreate            : %d (%s)\n", (int)rc, hipfftResultString(rc));
+    }
+  }
+
+  // 3. hipfftPlan1d (simplest API, batch = howmany)
+  {
+    hipfftHandle p;
+    hipfftResult rv = hipfftPlan1d(&p, G, HIPFFT_Z2Z, (int)howmany);
+    printf("  hipfftPlan1d     create : %d (%s)\n", (int)rv, hipfftResultString(rv));
+    if (rv == HIPFFT_SUCCESS) {
+      rv = hipfftExecZ2Z(p, dbuf, dbuf, HIPFFT_FORWARD);
+      hipDeviceSynchronize();
+      printf("  hipfftPlan1d     execFwd: %d (%s)\n", (int)rv, hipfftResultString(rv));
+      hipfftDestroy(p);
+    }
+  }
+
+  hipFree(dbuf);
+  printf("\n");
+}
+
+int main(void) {
+  // Print HIP device info
+  int device = 0;
+  hipGetDevice(&device);
+  hipDeviceProp_t prop;
+  hipGetDeviceProperties(&prop, device);
+  printf("Device %d: %s  warpSize=%d\n\n", device, prop.name, prop.warpSize);
+
+#ifdef hipfftVersionMinor
+  printf("hipFFT version: %d.%d.%d\n\n",
+         hipfftVersionMajor, hipfftVersionMinor, hipfftVersionPatch);
+#endif
+
+  // Original sweep with small howmany (these passed first time)
+  printf("=== Small howmany (original sweep) ===\n\n");
+  for (int G : {4, 8, 12, 16, 24, 32, 48, 64})
+    tryPlanAndExec(G, 512);
+
+  // Grid-realistic howmany values derived from actual lattice geometries.
+  // howmany = Ncomp * product(ldimensions[d] for d != dim)
+  // For LatticeComplexD: Ncomp=1.
+  printf("=== Grid-realistic parameters ===\n\n");
+
+  // --grid 16.16.16.16  4D FFT  (KNOWN TO FAIL in Grid)
+  // Each dim: G=16, Nperp=16^3=4096
+  tryPlanAndExec(16, 4096);
+
+  // --grid 32.32.32.32  4D FFT  (KNOWN TO SUCCEED in Grid)
+  // Each dim: G=32, Nperp=32^3=32768
+  tryPlanAndExec(32, 32768);
+
+  // --grid 32.32.32.32 Ls=8  5D DWF FFT  (KNOWN TO FAIL on dim 0 in Grid)
+  // dim 0: G=8,  Nperp=32^4=1048576
+  tryPlanAndExec(8, 1048576);
+  // dim 1-4: G=32, Nperp=8*32^3=262144
+  tryPlanAndExec(32, 262144);
+
+  // Extra intermediate cases to bracket the failure
+  tryPlanAndExec(16, 1024);
+  tryPlanAndExec(16, 2048);
+  tryPlanAndExec(16, 8192);
+  tryPlanAndExec(8,  4096);
+  tryPlanAndExec(8,  65536);
+  tryPlanAndExec(8,  262144);
+
+  return 0;
+}

tests/debug/Test_reduction.cc

@@ -73,36 +73,7 @@ void testReduction(GridCartesian *grid, GridParallelRNG &rng,
   Field field(grid);
 
   //--------------------------------------------------------------------
-  // 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.
+  // a) Timing: Niter timed calls reporting us/call and GB/s.
   //--------------------------------------------------------------------
 #if defined(GRID_CUDA) || defined(GRID_HIP) || defined(GRID_SYCL)
   {
@@ -113,38 +84,28 @@ 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_old(&v[0], osites);
+      for (int i = 0; i < Nwarm; i++) sum_gpu(&v[0], osites);
     }
 
-    RealD t_new, t_old;
+    RealD t_new;
     {
       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_new = bytes / (t_new / Niter) * 1e-3;
-    RealD GBs_old = bytes / (t_old / Niter) * 1e-3;
+    RealD bytes = (RealD)osites * sizeof(vobj);
+    RealD GBs   = bytes / (t_new / Niter) * 1e-3;
 
     std::cout << GridLogMessage << name << " timing (" << Niter << " calls):" << std::endl;
     std::cout << GridLogMessage
-              << "  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;
+              << "  sum_gpu " << t_new/Niter << " us   " << GBs << " GB/s" << std::endl;
   }
 #endif
 
   //--------------------------------------------------------------------
-  // d) Constant field via field = 1.0.
+  // b) 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)