tests/debug: extend hipfft reproducer with Grid-realistic howmany and exec tests

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

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
 
@@ -28,6 +28,10 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 #ifndef _GRID_FFT_H_
 #define _GRID_FFT_H_
 
+#include <any>
+#include <functional>
+#include <typeindex>
+
 #ifdef GRID_CUDA
 #include <cufft.h>
 #endif
@@ -65,17 +69,22 @@ 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) {
+    // hipfftPlanMany (one-step) triggers HIPFFT_PARSE_ERROR (12) on some
+    // ROCm versions.  The two-step hipfftCreate + hipfftMakePlanMany is
+    // more robust across ROCm releases.
     FFTW_plan p;
-    auto rv = hipfftPlanMany(&p,rank,n,n,istride,idist,n,ostride,odist,HIPFFT_Z2Z,howmany);
+    size_t workSize;
+    auto rc = hipfftCreate(&p);
+    GRID_ASSERT(rc==HIPFFT_SUCCESS);
+    auto rv = hipfftMakePlanMany(p,rank,n,nullptr,istride,idist,nullptr,ostride,odist,HIPFFT_Z2Z,howmany,&workSize);
     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 +92,28 @@ 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);
+    size_t workSize;
+    auto rc = hipfftCreate(&p);
+    GRID_ASSERT(rc==HIPFFT_SUCCESS);
+    auto rv = hipfftMakePlanMany(p,rank,n,nullptr,istride,idist,nullptr,ostride,odist,HIPFFT_C2C,howmany,&workSize);
     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 +121,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 +132,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 +181,314 @@ 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;
+
+  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;}    
+  GridCartesian *_grid;
 
-  FFT ( GridCartesian * grid ) 
-  {
-    flops=0;
-    usec =0;
+  // Type-erased plan entry.  The handle is recovered via
+  // std::any_cast<FFTW<scalar>::FFTW_plan> inside FFT_dim, which knows the
+  // scalar type at compile time.
+  struct PlanEntry {
+    std::any              handle;
+    std::function<void()> destroy;
   };
-    
-  ~FFT ( void)  {
-    //    delete sgrid;
-  }
-    
-  template<class vobj>
-  void FFT_dim_mask(Lattice<vobj> &result,const Lattice<vobj> &source,Coordinate mask,int sign){
 
-    //    vgrid=result.Grid();
-    //    conformable(result.Grid(),vgrid);
-    //    conformable(source.Grid(),vgrid);
+  std::vector<PlanEntry> forward_plans;   // size Nd when populated, 0 otherwise
+  std::vector<PlanEntry> backward_plans;
+  std::type_index        _plan_type { typeid(void) }; // vobj type plans were built for
+
+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; }
+
+  FFT(GridCartesian *grid) : _grid(grid), flops(0), usec(0) {}
+
+  ~FFT() {
+    if (forward_plans.size() > 0) PlanDestroy();
+  }
+
+  // Explicitly pre-create and cache plans for all Nd dimensions.
+  // Optional: FFT_dim will call this lazily on first use if not called.
+  // Asserts that no plans already exist; call PlanDestroy first to re-create.
+  template<class vobj>
+  void PlanCreate() {
+    GRID_ASSERT(forward_plans.size() == 0);
+
+    typedef typename vobj::scalar_type   scalar;
+    typedef typename vobj::scalar_object sobj;
+    typedef typename FFTW<scalar>::FFTW_scalar FFTW_scalar;
+    typedef typename FFTW<scalar>::FFTW_plan   FFTW_plan;
+
+    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};
+
+      // GPU backends (cuFFT/hipFFT) ignore the buffer pointer at plan creation.
+      // CPU FFTW with FFTW_ESTIMATE inspects only alignment and never touches data.
+      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, FFTW_FORWARD, FFTW_ESTIMATE);
+        forward_plans[d] = { p, [p](){ FFTW<scalar>::fftw_destroy_plan(p); } };
+      }
+      {
+        FFTW_plan p = FFTW<scalar>::fftw_plan_many_dft(
+            1, n, howmany, buf, n, 1, G, buf, n, 1, G, FFTW_BACKWARD, FFTW_ESTIMATE);
+        backward_plans[d] = { p, [p](){ FFTW<scalar>::fftw_destroy_plan(p); } };
+      }
+    }
+
+    _plan_type = std::type_index(typeid(vobj));
+  }
+
+  void PlanDestroy() {
+    for (auto &e : forward_plans)  e.destroy();
+    for (auto &e : backward_plans) e.destroy();
+    forward_plans.resize(0);
+    backward_plans.resize(0);
+    _plan_type = std::type_index(typeid(void));
+  }
+
+  template<class vobj>
+  void FFT_dim_mask(Lattice<vobj> &result, const Lattice<vobj> &source, Coordinate mask, int sign) {
     const int Ndim = source.Grid()->Nd();
     Lattice<vobj> tmp = source;
-    for(int d=0;d<Ndim;d++){
-      if( mask[d] ) {
-	FFT_dim(result,tmp,d,sign);
-	tmp=result;
+    for (int d = 0; d < Ndim; d++) {
+      if (mask[d]) {
+        FFT_dim(result, tmp, d, sign);
+        tmp = result;
       }
     }
   }
 
   template<class vobj>
-  void FFT_all_dim(Lattice<vobj> &result,const Lattice<vobj> &source,int sign){
+  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);
+    Coordinate mask(Ndim, 1);
+    FFT_dim_mask(result, source, mask, sign);
   }
 
-
   template<class vobj>
-  void FFT_dim(Lattice<vobj> &result,const Lattice<vobj> &source,int dim, int sign){
+  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());
+    conformable(result.Grid(), source.Grid());
 
     int L = grid->_ldimensions[dim];
     int G = grid->_fdimensions[dim];
-      
-    Coordinate layout(Ndim,1);
-    
-    // Construct pencils
+
     typedef typename vobj::scalar_object sobj;
-    typedef typename vobj::scalar_type   scalar;
     typedef typename vobj::scalar_type   scalar_type;
     typedef typename vobj::vector_type   vector_type;
-      
-    //std::cout << "CPU view" << std::endl;
-    
-    typedef typename FFTW<scalar>::FFTW_scalar FFTW_scalar;
-    typedef typename FFTW<scalar>::FFTW_plan   FFTW_plan;
-      
-    int Ncomp = sizeof(sobj)/sizeof(scalar);
-    int64_t Nlow  = 1;
-    int64_t Nhigh = 1;
 
-    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 */
+    typedef typename FFTW<scalar_type>::FFTW_scalar FFTW_scalar;
+    typedef typename FFTW<scalar_type>::FFTW_plan   FFTW_plan;
+
+    int     Ncomp  = sizeof(sobj) / sizeof(scalar_type);
+    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_type> pgbuf(Nperp * Ncomp * G); // [perp][component][dim]
+    scalar_type *pgbuf_v = &pgbuf[0];
+
+    int rank    = 1;
+    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 idist = G, odist = G, istride = 1, ostride = 1;
     int *inembed = n, *onembed = n;
-      
-    scalar div;
-    if ( sign == backward ) div = 1.0/G;
-    else if ( sign == forward ) div = 1.0;
+
+    scalar_type 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();
+    // Populate cache on first call; subsequent calls check type consistency.
+    if (forward_plans.size() == 0) PlanCreate<vobj>();
+    GRID_ASSERT(forward_plans.size() == (size_t)Ndim);
+    GRID_ASSERT(std::type_index(typeid(vobj)) == _plan_type);
+
+    auto &plans = (sign == forward) ? forward_plans : backward_plans;
+    FFTW_plan p = std::any_cast<FFTW_plan>(plans[dim].handle);
+
+    double t_pencil = 0;
+    double t_fft    = 0;
+    double t_copy   = 0;
+    double t_shift  = 0;
+    double t_total  = -usecond();
+
+    // Barrel-shift gather: accumulate global pencil into pgbuf
     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 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];
-    
+    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);
+    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); // buffer lane
+        {
+          int lane = acceleratorSIMTlane(Nsimd);
 #else
-      for(int lane=0;lane<Nsimd;lane++) {
+        for (int lane = 0; lane < Nsimd; lane++) {
 #endif
-	Coordinate icoor;
-	Coordinate ocoor;
-	Coordinate pgcoor;
+          Coordinate icoor, ocoor, pgcoor;
+          Lexicographic::CoorFromIndex(icoor, lane, sdims);
+          Lexicographic::CoorFromIndex(ocoor, idx,  rdims);
 
-	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);
 
-	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;
-	}
+          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();
 
-      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();
+      if (p_idx != processors[dim] - 1) {
+        Lattice<vobj> temp(grid);
+        t_shift -= usecond();
+        temp = Cshift(result, dim, L); result = temp;
+        t_shift += usecond();
       }
     }
     t_pencil += usecond();
-      
-    FFTW_scalar *in = (FFTW_scalar *)pgbuf_v;
-    FFTW_scalar *out= (FFTW_scalar *)pgbuf_v;
+
+    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);
+    FFTW<scalar_type>::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;
+
+    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,{
+      autoView(r_v, result, AcceleratorWrite);
+      accelerator_for(idx, grid->oSites(), Nsimd, {
 #ifdef GRID_SIMT
-      {
-	int lane=acceleratorSIMTlane(Nsimd); // buffer lane
+        {
+          int lane = acceleratorSIMTlane(Nsimd);
 #else
-      for(int lane=0;lane<Nsimd;lane++) {
+        for (int lane = 0; lane < Nsimd; lane++) {
 #endif
-	Coordinate icoor(Ndim);
-	Coordinate ocoor(Ndim);
-	Coordinate pgcoor(Ndim);
+          Coordinate icoor(Ndim), ocoor(Ndim), pgcoor(Ndim);
+          Lexicographic::CoorFromIndex(icoor, lane, sdims);
+          Lexicographic::CoorFromIndex(ocoor, idx,  rdims);
 
-	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);
 
-	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);
-	}
-	
+          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;
+    result   = result * div;
+    t_insert += usecond();
+    t_total  += usecond();
 
-    t_insert +=usecond();
-    
-    // destroying plan
-    FFTW<scalar>::fftw_destroy_plan(p);
-
-    t_total +=usecond();
-
-    std::cout <<GridLogPerformance<< " FFT took   "<<t_total/1.0e6 <<" s" << std::endl;
-    std::cout <<GridLogPerformance<< " FFT pencil "<<t_pencil/1.0e6 <<" s" << std::endl;
-    std::cout <<GridLogPerformance<< "  of which copy "<<t_copy/1.0e6 <<" s" << std::endl;
-    std::cout <<GridLogPerformance<< "  of which shift"<<t_shift/1.0e6 <<" s" << std::endl;
-    std::cout <<GridLogPerformance<< " FFT kernels "<<t_fft/1.0e6 <<" s" << std::endl;
-    std::cout <<GridLogPerformance<< " FFT insert  "<<t_insert/1.0e6 <<" s" << std::endl;
-    
+    std::cout << GridLogPerformance << " FFT took    " << t_total/1.0e6  << " s" << std::endl;
+    std::cout << GridLogPerformance << " FFT pencil  " << t_pencil/1.0e6 << " s" << std::endl;
+    std::cout << GridLogPerformance << "  of which copy  " << t_copy/1.0e6  << " s" << std::endl;
+    std::cout << GridLogPerformance << "  of which shift " << t_shift/1.0e6 << " s" << std::endl;
+    std::cout << GridLogPerformance << " FFT kernels " << t_fft/1.0e6    << " s" << std::endl;
+    std::cout << GridLogPerformance << " FFT insert  " << t_insert/1.0e6 << " s" << std::endl;
   }
 };
 

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

skills/gpu-memory-performance.md

@@ -9,13 +9,26 @@ allowed-tools:
 
 # 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 **2**) 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.
+`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 = 2;   // <-- default
+uint32_t accelerator_threads = 8;
 ```
 
 With `accelerator_for(ss, osites, nsimd, ...)`, the launch is:
@@ -25,15 +38,27 @@ dim3 threads(nsimd, acceleratorThreads(), 1)
 dim3 blocks ((osites + acceleratorThreads() - 1) / acceleratorThreads(), 1, 1)
 ```
 
-For `nsimd=1` and the default `acceleratorThreads()=2`:
-- **2 threads per block** on a 64-thread AMD wavefront → **3% occupancy**
-- Expected bandwidth ≈ peak × 3% ≈ 50 GB/s on MI250X
+Total threads per block = `Nsimd × acceleratorThreads()`.  With `GEN_SIMD_WIDTH=64B` and `Nsimd=8` (fp32 / ComplexF):
 
-**Diagnostic**: observed bandwidth << peak, kernel time >> expected from data volume.  Check with `--accelerator-threads 16` or `--accelerator-threads 32` at runtime.  A large speedup confirms occupancy starvation.
+| `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 16` or 32 at the application level.  Note this affects every `accelerator_for` site in the binary.
+2. Temporarily acceptable: set `--accelerator-threads 32` at the application level.  Note this affects every `accelerator_for` site in the binary.
 
 ## LambdaApply Thread Mapping
 
@@ -49,7 +74,7 @@ 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.
 
-Consequence: for coalesced access from a `vobj` array (AoS layout, stride = `sizeof(vobj)` between adjacent sites), adjacent threads in a wavefront address the **same** site at different lanes, not adjacent sites.  With `Nsimd=1` (GPU scalar build), `threadIdx.x` is always 0 and provides no coalescing benefit at all.
+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
 
@@ -57,7 +82,7 @@ These are Grid's canonical way to read/write one SIMD lane from a vector type in
 
 ```cpp
 // accelerator_for(ss, osites, Nsimd, {
-//   lane = acceleratorSIMTlane(Nsimd) = threadIdx.x
+//   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)
 ```
@@ -66,8 +91,6 @@ For `vobj` aggregate types, `coalescedRead` calls `extractLane(lane, vobj)` whic
 
 For `vsimd` (raw SIMD vector) types, it casts to `scalar_type*` and indexes with `lane`.
 
-**When Nsimd=1** (GPU scalar build): `lane=0` always, so `coalescedRead`/`coalescedWrite` are effectively no-ops (direct read/write).  Coalescing must be achieved through the iteration structure instead.
-
 ## 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:
@@ -108,12 +131,12 @@ Pattern: use `getNumBlocksAndThreads` to pick `numThreads` and `numBlocks`:
 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 → near-100% wavefront occupancy
+// gives 64–256 threads/block → correct occupancy independent of acceleratorThreads()
 Integer smemSize = numThreads * sizeof(sobj);
 myKernel<<<numBlocks, numThreads, smemSize, computeStream>>>(args...);
 ```
 
-This gives 64–256 threads/block regardless of `acceleratorThreads()`.  Grid's `reduceKernel` uses this pattern and achieves ~400 GB/s on MI250X.
+Grid's `reduceKernel` uses this pattern and achieves ~400 GB/s on MI250X.
 
 ## Fused vs Staged HBM Access
 
@@ -161,21 +184,22 @@ __device__ void packReduceBlocks(
 }
 ```
 
-Launched with `getNumBlocksAndThreads` → 128–256 threads/block → correct occupancy without depending on `acceleratorThreads()`.
+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=1)
+## 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, staged | 10,080 | 470 | 126,909 | 50 |
-| acceleratorThreads=16, staged | 342 | 310 | 8,251 | 297 |
-| acceleratorThreads=16, fused | — | 349 | 4,584 | 546 |
+| 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 threads per block: `accelerator_for(ss, N, 1, ...)` with default `acceleratorThreads()=2` = 2 threads/block = 3% occupancy on AMD.  Try `--accelerator-threads 16` at runtime; if it helps a lot, occupancy is the problem.
-2. Check for bulk struct accumulation in registers (`Bundle b; for(...) b._internal[k] = ...;`).  Replace with per-element writes via `coalescedWrite`.
-3. Check for staged HBM access (pack → buffer → reduce).  Count the passes; fuse if ≥ 2 passes over the same data.
-4. For reduction kernels, always use `getNumBlocksAndThreads` rather than `accelerator_for` so thread count is independent of `acceleratorThreads()`.
+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()`.

tests/core/Test_fft.cc

@@ -113,6 +113,7 @@ int main (int argc, char ** argv)
   Cref= Cref - C;
   std::cout << " invertible check " << norm2(Cref)<<std::endl;
 
+  theFFT.PlanDestroy();
   Stilde=S;
   std::cout<<" Benchmarking FFT of LatticeSpinMatrix  "<<std::endl;
   theFFT.FFT_dim(Stilde,Stilde,0,FFT::forward); std::cout << theFFT.MFlops()<<" mflops "<<std::endl;

tests/core/Test_fftf.cc

@@ -95,6 +95,7 @@ int main (int argc, char ** argv)
   C=C-Ctilde;
   std::cout << "diff scalar "<<norm2(C) << std::endl;
 
+  theFFT.PlanDestroy();
   Stilde = S;
   theFFT.FFT_dim(Stilde,Stilde,0,FFT::forward); std::cout << theFFT.MFlops()<< " "<<theFFT.USec() <<std::endl;
   theFFT.FFT_dim(Stilde,Stilde,1,FFT::forward); std::cout << theFFT.MFlops()<< " "<<theFFT.USec() <<std::endl;

tests/debug/Test_hipfft_minimal.cc

@@ -0,0 +1,158 @@
+/*
+ * 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) 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;
+}