Updated rocm 7 compile for ORNL

evils

Grid/algorithms/FFT.h

@@ -75,7 +75,6 @@ public:
     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,9 +82,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); }
 };
 template<> struct FFTW<ComplexF> {
 public:
@@ -93,7 +90,6 @@ public:
   static const int backward=FFTW_BACKWARD;
   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,
@@ -105,7 +101,6 @@ public:
     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,7 +119,6 @@ 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,
@@ -137,15 +129,12 @@ public:
     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:
@@ -153,7 +142,6 @@ public:
   static const int backward=FFTW_BACKWARD;
   typedef cufftComplex 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,
@@ -164,15 +152,12 @@ public:
     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
 
@@ -190,13 +175,10 @@ public:
 				      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:
@@ -210,25 +192,19 @@ public:
 				      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:
-    
+struct FFTbase {
   double        flops;
   double        flops_call;
   uint64_t      usec;
-    
-public:
+  GridCartesian *_grid;
 
   static const int forward  = FFTW_FORWARD;
   static const int backward = FFTW_BACKWARD;
@@ -237,168 +213,87 @@ public:
   double MFlops(void) { return flops / usec; }
   double USec(void)   { return (double)usec; }
 
-  FFT ( GridCartesian * grid ) 
-  {
-    flops=0;
-    usec =0;
+  FFTbase(GridCartesian *grid) : _grid(grid), flops(0), flops_call(0), usec(0) {}
 };
 
-  ~FFT ( void)  {
-    //    delete sgrid;
-  }
-    
+// 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>
-  void FFT_dim_mask(Lattice<vobj> &result,const Lattice<vobj> &source,Coordinate mask,int sign){
-
-    //    vgrid=result.Grid();
-    //    conformable(result.Grid(),vgrid);
-    //    conformable(source.Grid(),vgrid);
-    const int Ndim = source.Grid()->Nd();
-    Lattice<vobj> tmp = source;
-    for(int d=0;d<Ndim;d++){
-      if( mask[d] ) {
-	FFT_dim(result,tmp,d,sign);
-	tmp=result;
-      }
-    }
-  }
-
-  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);
-  }
-
-
-  template<class vobj>
-  void FFT_dim(Lattice<vobj> &result,const Lattice<vobj> &source,int dim, int sign){
-    const int Ndim = source.Grid()->Nd();
-    GridBase *grid = source.Grid();
-    conformable(result.Grid(),source.Grid());
-
-    int L = grid->_ldimensions[dim];
-    int G = grid->_fdimensions[dim];
-      
-    Coordinate layout(Ndim,1);
-    
-    // Construct pencils
-    typedef typename vobj::scalar_object sobj;
+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;
-      
-    //std::cout << "CPU view" << std::endl;
-    
   typedef typename FFTW<scalar>::FFTW_scalar FFTW_scalar;
-    typedef typename FFTW<scalar>::FFTW_plan   FFTW_plan;
 
+  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;
-    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 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; // Layout is [perp][component][dim]
-    pgbuf.resize(Nperp*Ncomp*G);
+  deviceVector<scalar> pgbuf(Nperp * Ncomp * G);
   scalar *pgbuf_v = &pgbuf[0];
-      
-    int rank = 1;  /* 1d transforms */
-    int n[] = {G}; /* 1d transforms of length G */
   int howmany = Ncomp * Nperp;
-    int odist,idist,istride,ostride;
-    idist   = odist   = G;            /* Distance between consecutive FT's */
-    istride = ostride = 1;            /* Distance between two elements in the same FT */
-    int *inembed = n, *onembed = n;
 
   scalar div;
-    if ( sign == backward ) div = 1.0/G;
-    else if ( sign == forward ) div = 1.0;
+  if      (sign == FFTW_BACKWARD) div = 1.0 / G;
+  else if (sign == FFTW_FORWARD)  div = 1.0;
   else GRID_ASSERT(0);
 
-    double t_pencil=0;
-    double t_fft   =0;
+  double t_pencil = 0, t_fft = 0, t_copy = 0, t_shift = 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;
+  const Coordinate processors = grid->_processors;
 
-    Coordinate processors = grid->_processors;
   Coordinate pgdims(Ndim);
   pgdims[0] = G;
-    for(int d=0, dd=1;d<Ndim;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_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++) {
 #endif
-	Coordinate icoor;
-	Coordinate ocoor;
-	Coordinate pgcoor;
-
+      Coordinate icoor, ocoor, 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
+      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++) {
-	  int64_t pg_idx = pgidx + w*pgvol;
         stmp = getlane(from[w], lane);
-	  pgbuf_v[pg_idx] = stmp;
+        pgbuf_v[pgidx + w*pgvol] = stmp;
       }
 #ifdef GRID_SIMT
     }
@@ -406,9 +301,8 @@ public:
     }
 #endif
     });
-
     t_copy += usecond();
-      if (p != processors[dim] - 1) {
+    if (p_idx != processors[dim] - 1) {
       Lattice<vobj> temp(grid);
       t_shift -= usecond();
       temp = Cshift(result, dim, L); result = temp;
@@ -423,49 +317,35 @@ public:
   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
+      int lane = acceleratorSIMTlane(Nsimd);
 #else
     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);
-
       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
+      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++) {
-	  int64_t pg_idx = pgidx + w*pgvol;
-	  stmp = pgbuf_v[pg_idx];
+        stmp = pgbuf_v[pgidx + w*pgvol];
         putlane(to[w], stmp, lane);
       }
-	
 #ifdef GRID_SIMT
     }
 #else
@@ -473,14 +353,8 @@ public:
 #endif
     });
   }
-
   result = result * div;
-
   t_insert += usecond();
-    
-    // destroying plan
-    FFTW<scalar>::fftw_destroy_plan(p);
-
   t_total  += usecond();
 
   std::cout << GridLogPerformance << " FFT took   " << t_total/1.0e6  << " s" << std::endl;
@@ -489,7 +363,130 @@ public:
   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;
+      }
+    }
+  }
+
+  template<class vobj>
+  void FFT_all_dim(Lattice<vobj> &result, const Lattice<vobj> &source, int sign) {
+    Coordinate mask(_grid->Nd(), 1);
+    FFT_dim_mask(result, source, mask, sign);
+  }
+
+  template<class vobj>
+  void FFT_dim(Lattice<vobj> &result, const Lattice<vobj> &source, int dim, int sign) {
+    GRID_ASSERT(source.Grid() == _grid);
+    GRID_ASSERT(result.Grid() == _grid);
+    conformable(result.Grid(), source.Grid());
+
+    typedef typename vobj::scalar_type   scalar;
+    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();
+    int G     = _grid->_fdimensions[dim];
+    int Ncomp = sizeof(sobj) / sizeof(scalar);
+    int64_t Nperp = 1;
+    for (int d = 0; d < Ndim; d++)
+      if (d != dim) Nperp *= _grid->_ldimensions[d];
+    int n[]     = {G};
+    int howmany = Ncomp * Nperp;
+
+    deviceVector<scalar> dummy(2);
+    FFTW_scalar *buf = (FFTW_scalar *)&dummy[0];
+    FFTW_plan p = FFTW<scalar>::fftw_plan_many_dft(1, n, howmany,
+                                                    buf, n, 1, G,
+                                                    buf, n, 1, G,
+                                                    sign, FFTW_ESTIMATE);
+    FFT_dim_execute(result, source, dim, sign, p, _grid, flops, flops_call, usec);
+    FFTW<scalar>::fftw_destroy_plan(p);
+  }
+};
+
+template<class vobj>
+class PlannedFFT : public FFTbase {
+private:
+  typedef typename vobj::scalar_type   scalar;
+  typedef typename vobj::scalar_object sobj;
+  typedef typename vobj::vector_type   vector_type;
+  typedef typename FFTW<scalar>::FFTW_scalar FFTW_scalar;
+  typedef typename FFTW<scalar>::FFTW_plan   FFTW_plan;
+
+  std::vector<FFTW_plan> forward_plans;
+  std::vector<FFTW_plan> backward_plans;
+
+  void PlanCreate() {
+    const int Ndim = _grid->Nd();
+    forward_plans.resize(Ndim);
+    backward_plans.resize(Ndim);
+
+    for (int d = 0; d < Ndim; d++) {
+      int G       = _grid->_fdimensions[d];
+      int Ncomp   = sizeof(sobj) / sizeof(scalar);
+      int64_t Nperp = 1;
+      for (int dd = 0; dd < Ndim; dd++)
+        if (dd != d) Nperp *= _grid->_ldimensions[dd];
+      int howmany = Ncomp * (int)Nperp;
+      int n[]     = {G};
+
+      deviceVector<scalar> dummy(2);
+      FFTW_scalar *buf = (FFTW_scalar *)&dummy[0];
+
+      forward_plans[d]  = FFTW<scalar>::fftw_plan_many_dft(1, n, howmany, buf, n, 1, G, buf, n, 1, G, FFTW_FORWARD,  FFTW_ESTIMATE);
+      backward_plans[d] = FFTW<scalar>::fftw_plan_many_dft(1, n, howmany, buf, n, 1, G, buf, n, 1, G, FFTW_BACKWARD, FFTW_ESTIMATE);
+    }
+  }
+
+  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/lattice/Lattice_reduction_gpu.h

@@ -198,7 +198,7 @@ __global__ void reduceKernel(const vobj *lat, sobj *buffer, Iterator n) {
 // Possibly promote to double and sum
 /////////////////////////////////////////////////////////////////////////////////////////////////////////
 
-#define GRID_REDUCTION_TIMING
+#undef GRID_REDUCTION_TIMING
 
 template <class vobj>
 inline typename vobj::scalar_objectD sumD_gpu_small(const vobj *lat, Integer osites)
@@ -230,7 +230,7 @@ inline typename vobj::scalar_objectD sumD_gpu_small(const vobj *lat, Integer osi
   acceleratorCopyFromDevice(buffer_v,&result,sizeof(result));
 #ifdef GRID_REDUCTION_TIMING
   t_d2h += usecond();
-  std::cout << GridLogMessage << "  sumD_gpu_small"
+  std::cout << GridLogDebug << "  sumD_gpu_small"
             << " sizeof(sobj)=" << sizeof(sobj)
             << " blocks=" << numBlocks << " threads=" << numThreads
             << " kernel+barrier=" << t_kernel << " us"
@@ -362,7 +362,7 @@ inline void sumD_gpu_reduce_words(const vobj *lat, Integer osites,
   acceleratorCopyFromDevice(buffer_v, &result, sizeof(result));
 #ifdef GRID_REDUCTION_TIMING
   t_d2h += usecond();
-  std::cout << GridLogMessage << " sumD_gpu_reduce_words R=" << R
+  std::cout << GridLogDebug << " sumD_gpu_reduce_words R=" << R
             << " base=" << base
             << " kernel=" << t_kernel << " D2H=" << t_d2h << " us" << std::endl;
 #endif
@@ -391,7 +391,7 @@ inline typename vobj::scalar_objectD sumD_gpu_large(const vobj *lat, Integer osi
   while (w      <  words) { sumD_gpu_reduce_words< 1>(lat, osites, ret_p, w); w +=  1; }
 #ifdef GRID_REDUCTION_TIMING
   t_large += usecond();
-  std::cout << GridLogMessage << "sumD_gpu_large"
+  std::cout << GridLogDebug << "sumD_gpu_large"
             << " sizeof(sobjD)=" << sizeof(sobjD)
             << " words=" << words << " total=" << t_large << " us" << std::endl;
 #endif

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

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

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

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:

systems/Frontier/config-command

@@ -13,8 +13,8 @@ CLIME=`spack find --paths c-lime@2-3-9 | grep c-lime| cut -c 15-`
 --with-mpfr=/opt/cray/pe/gcc/mpfr/3.1.4/ \
 --disable-fermion-reps \
 CXX=hipcc MPICXX=mpicxx \
-CXXFLAGS="-fPIC -I${ROCM_PATH}/include/ -I${MPICH_DIR}/include -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,10 @@
 
 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
 
-#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

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

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