Fixing inconsistency in configure.ac

worth maintaining

Grid/algorithms/FFT.h

@@ -28,10 +28,6 @@ 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
@@ -74,14 +70,8 @@ public:
 				      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;
-    size_t workSize;
+    auto rv = hipfftPlanMany(&p,rank,n,n,istride,idist,n,ostride,odist,HIPFFT_Z2Z,howmany);
-    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;
   }
@@ -107,10 +97,7 @@ public:
 				      int ostride, int odist,
 				      int sign, unsigned flags) {
     FFTW_plan p;
-    size_t workSize;
+    auto rv = hipfftPlanMany(&p,rank,n,n,istride,idist,n,ostride,odist,HIPFFT_C2C,howmany);
-    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;
   }
@@ -213,28 +200,12 @@ public:
 #endif
 #endif
 
-class FFT {
+struct FFTbase {
-private:
+  double        flops;
-
+  double        flops_call;
-  double   flops;
+  uint64_t      usec;
-  double   flops_call;
-  uint64_t usec;
   GridCartesian *_grid;
 
-  // 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;
-  };
-
-  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;
 
@@ -242,70 +213,166 @@ public:
   double MFlops(void) { return flops / usec; }
   double USec(void)   { return (double)usec; }
 
-  FFT(GridCartesian *grid) : _grid(grid), flops(0), usec(0) {}
+  FFTbase(GridCartesian *grid) : _grid(grid), flops(0), flops_call(0), usec(0) {}
+};
 
-  ~FFT() {
+// Barrel-shift gather, FFT execute, and insert. Called by both FFT and PlannedFFT.
-    if (forward_plans.size() > 0) PlanDestroy();
+// The caller is responsible for plan acquisition and destruction.
-  }
+template<class vobj>
+static void FFT_dim_execute(
+    Lattice<vobj>       &result,
+    const Lattice<vobj> &source,
+    int dim, int sign,
+    typename FFTW<typename vobj::scalar_type>::FFTW_plan p,
+    GridCartesian *grid,
+    double &flops, double &flops_call, uint64_t &usec)
+{
+  typedef typename vobj::scalar_type   scalar;
+  typedef typename vobj::scalar_object sobj;
+  typedef typename vobj::scalar_type   scalar_type;
+  typedef typename vobj::vector_type   vector_type;
+  typedef typename FFTW<scalar>::FFTW_scalar FFTW_scalar;
 
-  // Explicitly pre-create and cache plans for all Nd dimensions.
+  const int Ndim = grid->Nd();
-  // Optional: FFT_dim will call this lazily on first use if not called.
+  int L = grid->_ldimensions[dim];
-  // Asserts that no plans already exist; call PlanDestroy first to re-create.
+  int G = grid->_fdimensions[dim];
-  template<class vobj>
+  int Ncomp = sizeof(sobj) / sizeof(scalar);
-  void PlanCreate() {
+  int64_t Nlow = 1, Nhigh = 1;
-    GRID_ASSERT(forward_plans.size() == 0);
+  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;
 
-    typedef typename vobj::scalar_type   scalar;
+  deviceVector<scalar> pgbuf(Nperp * Ncomp * G);
-    typedef typename vobj::scalar_object sobj;
+  scalar *pgbuf_v = &pgbuf[0];
-    typedef typename FFTW<scalar>::FFTW_scalar FFTW_scalar;
+  int howmany = Ncomp * Nperp;
-    typedef typename FFTW<scalar>::FFTW_plan   FFTW_plan;
 
-    const int Ndim = _grid->Nd();
+  scalar div;
-    forward_plans.resize(Ndim);
+  if      (sign == FFTW_BACKWARD) div = 1.0 / G;
-    backward_plans.resize(Ndim);
+  else if (sign == FFTW_FORWARD)  div = 1.0;
+  else GRID_ASSERT(0);
 
-    for (int d = 0; d < Ndim; d++) {
+  double t_pencil = 0, t_fft = 0, t_copy = 0, t_shift = 0;
-      int G     = _grid->_fdimensions[d];
+  double t_total  = -usecond();
-      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.
+  result = source;
-      // CPU FFTW with FFTW_ESTIMATE inspects only alignment and never touches data.
+  int pc = grid->_processor_coor[dim];
-      // Use a host stack buffer: a device allocation here triggers a rocFFT RTC
-      // bug on ROCm 7 that causes plan creation to fail for small transform sizes.
-      scalar stack_dummy[2] = {};
-      FFTW_scalar *buf = (FFTW_scalar *)stack_dummy;
 
-      {
+  const Coordinate ldims      = grid->_ldimensions;
-        FFTW_plan p = FFTW<scalar>::fftw_plan_many_dft(
+  const Coordinate rdims      = grid->_rdimensions;
-            1, n, howmany, buf, n, 1, G, buf, n, 1, G, FFTW_FORWARD, FFTW_ESTIMATE);
+  const Coordinate sdims      = grid->_simd_layout;
-        forward_plans[d] = { p, [p](){ FFTW<scalar>::fftw_destroy_plan(p); } };
+  const Coordinate processors = grid->_processors;
-      }
+
-      {
+  Coordinate pgdims(Ndim);
-        FFTW_plan p = FFTW<scalar>::fftw_plan_many_dft(
+  pgdims[0] = G;
-            1, n, howmany, buf, n, 1, G, buf, n, 1, G, FFTW_BACKWARD, FFTW_ESTIMATE);
+  for (int d = 0, dd = 1; d < Ndim; d++)
-        backward_plans[d] = { p, [p](){ FFTW<scalar>::fftw_destroy_plan(p); } };
+    if (d != dim) pgdims[dd++] = ldims[d];
+  int64_t pgvol = 1;
+  for (int d = 0; d < Ndim; d++) pgvol *= pgdims[d];
+
+  const int Nsimd = vobj::Nsimd();
+  t_pencil = -usecond();
+  for (int p_idx = 0; p_idx < processors[dim]; p_idx++) {
+    t_copy -= usecond();
+    autoView(r_v, result, AcceleratorRead);
+    accelerator_for(idx, grid->oSites(), vobj::Nsimd(), {
+#ifdef GRID_SIMT
+    {
+      int lane = acceleratorSIMTlane(Nsimd);
+#else
+    for (int lane = 0; lane < Nsimd; lane++) {
+#endif
+      Coordinate icoor, ocoor, pgcoor;
+      Lexicographic::CoorFromIndex(icoor, lane, sdims);
+      Lexicographic::CoorFromIndex(ocoor, idx,  rdims);
+      pgcoor[0] = ocoor[dim] + icoor[dim]*rdims[dim] + ((pc+p_idx)%processors[dim])*L;
+      for (int d = 0, dd = 1; d < Ndim; d++)
+        if (d != dim) { pgcoor[dd] = ocoor[d] + icoor[d]*rdims[d]; dd++; }
+      int64_t pgidx;
+      Lexicographic::IndexFromCoor(pgcoor, pgidx, pgdims);
+      vector_type *from = (vector_type *)&r_v[idx];
+      scalar_type  stmp;
+      for (int w = 0; w < Ncomp; w++) {
+        stmp = getlane(from[w], lane);
+        pgbuf_v[pgidx + w*pgvol] = stmp;
       }
+#ifdef GRID_SIMT
+    }
+#else
+    }
+#endif
+    });
+    t_copy += usecond();
+    if (p_idx != processors[dim] - 1) {
+      Lattice<vobj> temp(grid);
+      t_shift -= usecond();
+      temp = Cshift(result, dim, L); result = temp;
+      t_shift += usecond();
     }
-
-    _plan_type = std::type_index(typeid(vobj));
   }
+  t_pencil += usecond();
 
-  void PlanDestroy() {
+  FFTW_scalar *in  = (FFTW_scalar *)pgbuf_v;
-    for (auto &e : forward_plans)  e.destroy();
+  FFTW_scalar *out = (FFTW_scalar *)pgbuf_v;
-    for (auto &e : backward_plans) e.destroy();
+  t_fft = -usecond();
-    forward_plans.resize(0);
+  FFTW<scalar>::fftw_execute_dft(p, in, out, sign);
-    backward_plans.resize(0);
+  t_fft += usecond();
-    _plan_type = std::type_index(typeid(void));
+
+  flops_call = 5.0 * howmany * G * log2(G);
+  usec  = t_fft;
+  flops = flops_call;
+
+  result = Zero();
+  double t_insert = -usecond();
+  {
+    autoView(r_v, result, AcceleratorWrite);
+    accelerator_for(idx, grid->oSites(), Nsimd, {
+#ifdef GRID_SIMT
+    {
+      int lane = acceleratorSIMTlane(Nsimd);
+#else
+    for (int lane = 0; lane < Nsimd; lane++) {
+#endif
+      Coordinate icoor(Ndim), ocoor(Ndim), pgcoor(Ndim);
+      Lexicographic::CoorFromIndex(icoor, lane, sdims);
+      Lexicographic::CoorFromIndex(ocoor, idx,  rdims);
+      pgcoor[0] = ocoor[dim] + icoor[dim]*rdims[dim] + pc*L;
+      for (int d = 0, dd = 1; d < Ndim; d++)
+        if (d != dim) { pgcoor[dd] = ocoor[d] + icoor[d]*rdims[d]; dd++; }
+      int64_t pgidx;
+      Lexicographic::IndexFromCoor(pgcoor, pgidx, pgdims);
+      vector_type *to  = (vector_type *)&r_v[idx];
+      scalar_type  stmp;
+      for (int w = 0; w < Ncomp; w++) {
+        stmp = pgbuf_v[pgidx + w*pgvol];
+        putlane(to[w], stmp, lane);
+      }
+#ifdef GRID_SIMT
+    }
+#else
+    }
+#endif
+    });
   }
+  result = result * div;
+  t_insert += usecond();
+  t_total  += usecond();
+
+  std::cout << GridLogPerformance << " FFT took   " << t_total/1.0e6  << " s" << std::endl;
+  std::cout << GridLogPerformance << " FFT pencil " << t_pencil/1.0e6 << " s" << std::endl;
+  std::cout << GridLogPerformance << "  of which copy " << t_copy/1.0e6  << " s" << std::endl;
+  std::cout << GridLogPerformance << "  of which shift" << t_shift/1.0e6 << " s" << std::endl;
+  std::cout << GridLogPerformance << " FFT kernels " << t_fft/1.0e6    << " s" << std::endl;
+  std::cout << GridLogPerformance << " FFT insert  " << t_insert/1.0e6 << " s" << std::endl;
+}
+
+class FFT : public FFTbase {
+public:
+  FFT(GridCartesian *grid) : FFTbase(grid) {}
+  ~FFT() {}
 
   template<class vobj>
   void FFT_dim_mask(Lattice<vobj> &result, const Lattice<vobj> &source, Coordinate mask, int sign) {
-    const int Ndim = source.Grid()->Nd();
+    const int Ndim = _grid->Nd();
     Lattice<vobj> tmp = source;
     for (int d = 0; d < Ndim; d++) {
       if (mask[d]) {
@@ -317,180 +384,109 @@ public:
 
   template<class vobj>
   void FFT_all_dim(Lattice<vobj> &result, const Lattice<vobj> &source, int sign) {
-    const int Ndim = source.Grid()->Nd();
+    Coordinate mask(_grid->Nd(), 1);
-    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();
+    GRID_ASSERT(source.Grid() == _grid);
-    GridBase *grid = source.Grid();
+    GRID_ASSERT(result.Grid() == _grid);
     conformable(result.Grid(), source.Grid());
 
-    int L = grid->_ldimensions[dim];
+    typedef typename vobj::scalar_type   scalar;
-    int G = grid->_fdimensions[dim];
-
     typedef typename vobj::scalar_object sobj;
-    typedef typename vobj::scalar_type   scalar_type;
+    typedef typename FFTW<scalar>::FFTW_scalar FFTW_scalar;
-    typedef typename vobj::vector_type   vector_type;
+    typedef typename FFTW<scalar>::FFTW_plan   FFTW_plan;
 
-    typedef typename FFTW<scalar_type>::FFTW_scalar FFTW_scalar;
+    const int Ndim = _grid->Nd();
-    typedef typename FFTW<scalar_type>::FFTW_plan   FFTW_plan;
+    int G     = _grid->_fdimensions[dim];
-
+    int Ncomp = sizeof(sobj) / sizeof(scalar);
-    int     Ncomp  = sizeof(sobj) / sizeof(scalar_type);
+    int64_t Nperp = 1;
-    int64_t Nlow   = 1;
+    for (int d = 0; d < Ndim; d++)
-    int64_t Nhigh  = 1;
+      if (d != dim) Nperp *= _grid->_ldimensions[d];
-    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 idist = G, odist = G, istride = 1, ostride = 1;
-    int *inembed = n, *onembed = n;
 
-    scalar_type div;
+    deviceVector<scalar> dummy(2);
-    if      (sign == backward) div = 1.0 / G;
+    FFTW_scalar *buf = (FFTW_scalar *)&dummy[0];
-    else if (sign == forward)  div = 1.0;
+    FFTW_plan p = FFTW<scalar>::fftw_plan_many_dft(1, n, howmany,
-    else GRID_ASSERT(0);
+                                                    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);
+  }
+};
 
-    // Populate cache on first call; subsequent calls check type consistency.
+template<class vobj>
-    if (forward_plans.size() == 0) PlanCreate<vobj>();
+class PlannedFFT : public FFTbase {
-    GRID_ASSERT(forward_plans.size() == (size_t)Ndim);
+private:
-    GRID_ASSERT(std::type_index(typeid(vobj)) == _plan_type);
+  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;
 
-    auto &plans = (sign == forward) ? forward_plans : backward_plans;
+  std::vector<FFTW_plan> forward_plans;
-    FFTW_plan p = std::any_cast<FFTW_plan>(plans[dim].handle);
+  std::vector<FFTW_plan> backward_plans;
 
-    double t_pencil = 0;
+  void PlanCreate() {
-    double t_fft    = 0;
+    const int Ndim = _grid->Nd();
-    double t_copy   = 0;
+    forward_plans.resize(Ndim);
-    double t_shift  = 0;
+    backward_plans.resize(Ndim);
-    double t_total  = -usecond();
 
-    // Barrel-shift gather: accumulate global pencil into pgbuf
+    for (int d = 0; d < Ndim; d++) {
-    result = source;
+      int G       = _grid->_fdimensions[d];
-    int pc = grid->_processor_coor[dim];
+      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};
 
-    const Coordinate ldims = grid->_ldimensions;
+      deviceVector<scalar> dummy(2);
-    const Coordinate rdims = grid->_rdimensions;
+      FFTW_scalar *buf = (FFTW_scalar *)&dummy[0];
-    const Coordinate sdims = grid->_simd_layout;
-    Coordinate processors  = grid->_processors;
 
-    Coordinate pgdims(Ndim);
+      forward_plans[d]  = FFTW<scalar>::fftw_plan_many_dft(1, n, howmany, buf, n, 1, G, buf, n, 1, G, FFTW_FORWARD,  FFTW_ESTIMATE);
-    pgdims[0] = G;
+      backward_plans[d] = FFTW<scalar>::fftw_plan_many_dft(1, n, howmany, buf, n, 1, G, buf, n, 1, G, FFTW_BACKWARD, FFTW_ESTIMATE);
-    for (int d = 0, 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();
+  void PlanDestroy() {
-    t_pencil = -usecond();
+    for (auto p : forward_plans)  FFTW<scalar>::fftw_destroy_plan(p);
-    for (int p_idx = 0; p_idx < processors[dim]; p_idx++) {
+    for (auto p : backward_plans) FFTW<scalar>::fftw_destroy_plan(p);
-      t_copy -= usecond();
+    forward_plans.clear();
-      autoView(r_v, result, AcceleratorRead);
+    backward_plans.clear();
-      accelerator_for(idx, grid->oSites(), vobj::Nsimd(), {
+  }
-#ifdef GRID_SIMT
-        {
-          int lane = acceleratorSIMTlane(Nsimd);
-#else
-        for (int lane = 0; lane < Nsimd; lane++) {
-#endif
-          Coordinate icoor, ocoor, pgcoor;
-          Lexicographic::CoorFromIndex(icoor, lane, sdims);
-          Lexicographic::CoorFromIndex(ocoor, idx,  rdims);
 
-          pgcoor[0] = ocoor[dim] + icoor[dim]*rdims[dim] + ((pc+p_idx)%processors[dim])*L;
+public:
-          for (int d = 0, dd = 1; d < Ndim; d++) {
+  PlannedFFT(GridCartesian *grid) : FFTbase(grid) { PlanCreate(); }
-            if (d != dim) { pgcoor[dd] = ocoor[d] + icoor[d]*rdims[d]; dd++; }
+  ~PlannedFFT() { PlanDestroy(); }
-          }
-          int64_t pgidx;
-          Lexicographic::IndexFromCoor(pgcoor, pgidx, pgdims);
 
-          vector_type   *from = (vector_type *)&r_v[idx];
+  void FFT_dim_mask(Lattice<vobj> &result, const Lattice<vobj> &source, Coordinate mask, int sign) {
-          scalar_type    stmp;
+    const int Ndim = _grid->Nd();
-          for (int w = 0; w < Ncomp; w++) {
+    Lattice<vobj> tmp = source;
-            stmp = getlane(from[w], lane);
+    for (int d = 0; d < Ndim; d++) {
-            pgbuf_v[pgidx + w*pgvol] = stmp;
+      if (mask[d]) {
-          }
+        FFT_dim(result, tmp, d, sign);
-#ifdef GRID_SIMT
+        tmp = result;
-        }
-#else
-        }
-#endif
-      });
-      t_copy += usecond();
-
-      if (p_idx != processors[dim] - 1) {
-        Lattice<vobj> temp(grid);
-        t_shift -= usecond();
-        temp = Cshift(result, dim, L); result = temp;
-        t_shift += usecond();
       }
     }
-    t_pencil += usecond();
+  }
 
-    FFTW_scalar *in  = (FFTW_scalar *)pgbuf_v;
+  void FFT_all_dim(Lattice<vobj> &result, const Lattice<vobj> &source, int sign) {
-    FFTW_scalar *out = (FFTW_scalar *)pgbuf_v;
+    Coordinate mask(_grid->Nd(), 1);
-    t_fft = -usecond();
+    FFT_dim_mask(result, source, mask, sign);
-    FFTW<scalar_type>::fftw_execute_dft(p, in, out, sign);
+  }
-    t_fft += usecond();
 
-    flops_call = 5.0 * howmany * G * log2(G);
+  void FFT_dim(Lattice<vobj> &result, const Lattice<vobj> &source, int dim, int sign) {
-    usec  = t_fft;
+    GRID_ASSERT(source.Grid() == _grid);
-    flops = flops_call;
+    GRID_ASSERT(result.Grid() == _grid);
-
+    GRID_ASSERT((int)forward_plans.size() == _grid->Nd());
-    result = Zero();
+    conformable(result.Grid(), source.Grid());
-    double t_insert = -usecond();
+    FFTW_plan p = (sign == forward ? forward_plans : backward_plans)[dim];
-    {
+    FFT_dim_execute(result, source, dim, sign, p, _grid, flops, flops_call, usec);
-      autoView(r_v, result, AcceleratorWrite);
-      accelerator_for(idx, grid->oSites(), Nsimd, {
-#ifdef GRID_SIMT
-        {
-          int lane = acceleratorSIMTlane(Nsimd);
-#else
-        for (int lane = 0; lane < Nsimd; lane++) {
-#endif
-          Coordinate icoor(Ndim), ocoor(Ndim), pgcoor(Ndim);
-          Lexicographic::CoorFromIndex(icoor, lane, sdims);
-          Lexicographic::CoorFromIndex(ocoor, idx,  rdims);
-
-          pgcoor[0] = ocoor[dim] + icoor[dim]*rdims[dim] + pc*L;
-          for (int d = 0, dd = 1; d < Ndim; d++) {
-            if (d != dim) { pgcoor[dd] = ocoor[d] + icoor[d]*rdims[d]; dd++; }
-          }
-          int64_t pgidx;
-          Lexicographic::IndexFromCoor(pgcoor, pgidx, pgdims);
-
-          vector_type *to   = (vector_type *)&r_v[idx];
-          scalar_type  stmp;
-          for (int w = 0; w < Ncomp; w++) {
-            stmp = pgbuf_v[pgidx + w*pgvol];
-            putlane(to[w], stmp, lane);
-          }
-#ifdef GRID_SIMT
-        }
-#else
-        }
-#endif
-      });
-    }
-
-    result   = result * div;
-    t_insert += usecond();
-    t_total  += usecond();
-
-    std::cout << GridLogPerformance << " FFT took    " << t_total/1.0e6  << " s" << std::endl;
-    std::cout << GridLogPerformance << " FFT pencil  " << t_pencil/1.0e6 << " s" << std::endl;
-    std::cout << GridLogPerformance << "  of which copy  " << t_copy/1.0e6  << " s" << std::endl;
-    std::cout << GridLogPerformance << "  of which shift " << t_shift/1.0e6 << " s" << std::endl;
-    std::cout << GridLogPerformance << " FFT kernels " << t_fft/1.0e6    << " s" << std::endl;
-    std::cout << GridLogPerformance << " FFT insert  " << t_insert/1.0e6 << " s" << std::endl;
   }
 };
 

Grid/allocator/MemoryManager.cc

@@ -63,12 +63,10 @@ void MemoryManager::PrintBytes(void)
   std::cout << " MemoryManager : "<<(total_device>>20)<<" accelerator Mbytes "<<std::endl;
   std::cout << " MemoryManager : "<<(total_host>>20)  <<" cpu         Mbytes "<<std::endl;
   uint64_t cacheBytes;
-  cacheBytes = CacheBytes[Cpu];
+  cacheBytes = HostCacheBytes();
-  std::cout << " MemoryManager : "<<(cacheBytes>>20) <<" cpu cache Mbytes "<<std::endl;
+  std::cout << " MemoryManager : "<<(cacheBytes>>20) <<" cpu alloc cache Mbytes "<<std::endl;
-  cacheBytes = CacheBytes[Acc];
+  cacheBytes = DeviceCacheBytes();
-  std::cout << " MemoryManager : "<<(cacheBytes>>20) <<" acc cache Mbytes "<<std::endl;
+  std::cout << " MemoryManager : "<<(cacheBytes>>20) <<" acc alloc cache Mbytes "<<std::endl;
-  cacheBytes = CacheBytes[Shared];
-  std::cout << " MemoryManager : "<<(cacheBytes>>20) <<" shared cache Mbytes "<<std::endl;
   
 #ifdef GRID_CUDA
   cuda_mem();

Grid/lattice/Lattice_base.h

@@ -289,7 +289,7 @@ public:
   ///////////////////////////////////////////
   // move constructor
   ///////////////////////////////////////////
-  Lattice(Lattice && r){ 
+  Lattice(Lattice && r) noexcept { 
     this->_grid = r.Grid();
     this->_odata      = r._odata;
     this->_odata_size = r._odata_size;
@@ -330,7 +330,7 @@ public:
   ///////////////////////////////////////////
   // Move assignment possible if same type
   ///////////////////////////////////////////
-  inline Lattice<vobj> & operator = (Lattice<vobj> && r){
+  inline Lattice<vobj> & operator = (Lattice<vobj> && r) noexcept {
 
     resize(0); // deletes if appropriate
     this->_grid       = r.Grid();

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
@@ -438,5 +438,11 @@ inline typename vobj::scalar_object sum_gpu_large(const vobj *lat, Integer osite
   result = sumD_gpu_large(lat,osites);
   return result;
 }
+template<class Word> Word checksum_gpu(Word *vec,uint64_t L)
+{
+  Word w;
+  bzero(&w,sizeof(w));
+  return w;
+}
 
 NAMESPACE_END(Grid);

Grid/lattice/Lattice_slicesum_core.h

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

Grid/perfmon/PerfCount.h

@@ -51,8 +51,8 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #endif
 #ifdef __x86_64__
 #ifdef GRID_CUDA
-//accelerator_inline uint64_t __rdtsc(void) {  return 0; }
+accelerator_inline uint64_t __rdtsc(void) {  return 0; }
-//accelerator_inline uint64_t __rdpmc(int ) {  return 0; }
+accelerator_inline uint64_t __rdpmc(int ) {  return 0; }
 #else
 #include <x86intrin.h>
 #endif
@@ -93,8 +93,7 @@ inline uint64_t cyclecount(void){
 }
 #elif defined __x86_64__
 inline uint64_t cyclecount(void){ 
-  uint64_t ret = __rdtsc();
+  return (uint64_t)0;
-  return (uint64_t)ret;
 }
 #else
 

Grid/serialisation/XmlIO.h

@@ -184,11 +184,29 @@ namespace Grid
       std::cout << GridLogWarning << "XML: cannot open node '" << s << "'";
       std::cout << std::endl;
     } else {
-      for(unsigned int i = 0; node_.child("elem"); )
+      // Consumed children are *removed* rather than renamed to "elem-done".
+      //
+      // pugixml's node_.child("elem") scans the child list from the start.
+      // Leaving consumed nodes in place under a new name meant iteration i
+      // had to skip i dead siblings, and the lookup happens three times per
+      // element (loop condition, read, rename) -- so reading an N-element
+      // vector cost O(N^2) node visits.  Measured on this reader: 11 s for
+      // 32k elements, 278 s for 131k, i.e. hours for the ~500k elements a
+      // 512x512 complex matrix produces.  Removing the node instead keeps
+      // the next "elem" as the first child, making each lookup O(1).
+      //
+      // Counting first also lets the vector be sized once instead of grown
+      // one element at a time.
+      size_t n = 0;
+      for (pugi::xml_node e = node_.child("elem"); e; e = e.next_sibling("elem"))
       {
-        output.resize(i + 1);
+        n++;
-        read("elem", output[i++]);
+      }
-        node_.child("elem").set_name("elem-done");
+      output.resize(n);
+      for (size_t i = 0; i < n; ++i)
+      {
+        read("elem", output[i]);
+        node_.remove_child(node_.child("elem"));
       }
       pop();
     }
@@ -225,7 +243,11 @@ namespace Grid
       {
         NumElements--;
         read("elem", x);
-        node_.child("elem").set_name("elem-done");
+        // Remove rather than rename, for the same O(N^2) reason as the
+        // std::vector reader above.  Only "elem" nodes are removed; the
+        // tensor-level renaming below is left alone because push(sName)
+        // relies on it to find the next unvisited tensor block.
+        node_.remove_child(node_.child("elem"));
         int i;
         for( i = Rank - 1 ; i != -1 && ++MyIndex[i] == dim[i] ; i-- )
           MyIndex[i] = 0;

Grid/simd/Simd.h

@@ -113,6 +113,14 @@ accelerator_inline RealD    adj(const RealD  & r){ return r; }
 accelerator_inline ComplexD adj(const ComplexD& r){ return(conjugate(r)); }
 accelerator_inline ComplexF adj(const ComplexF& r ){ return(conjugate(r)); }
 
+#if defined(GRID_CUDA) || defined(GRID_HIP)
+//Provide for convenience
+inline std::complex<double> conjugate(const std::complex<double>& r){ return(conj(r)); }
+inline std::complex<float>  conjugate(const std::complex<float>& r) { return(conj(r)); }
+inline std::complex<double> adj(const std::complex<double>& r)      { return(conj(r)); }
+inline std::complex<float>  adj(const std::complex<float>& r)       { return(conj(r)); }
+#endif
+
 accelerator_inline RealF real(const RealF  & r){ return r; }
 accelerator_inline RealD real(const RealD  & r){ return r; }
 accelerator_inline RealF real(const ComplexF  & r){ return r.real(); }

Grid/threads/Accelerator.h

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

configure.ac

@@ -174,7 +174,7 @@ esac
 ############### fermions
 AC_ARG_ENABLE([fermion-instantiations],
      [AS_HELP_STRING([--enable-fermion-instantiations=yes|no],[enable fermion instantiations])],
-     [ac_FERMION_REPS=${enable_fermion_instantiations}], [ac_FERMION_INSTANTIATIONS=yes])
+     [ac_FERMION_INSTANTIATIONS=${enable_fermion_instantiations}], [ac_FERMION_INSTANTIATIONS=yes])
 
 AM_CONDITIONAL(BUILD_FERMION_INSTANTIATIONS, [ test "${ac_FERMION_INSTANTIATIONS}X" == "yesX" ])
 

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/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/Booster/config-command

@@ -8,7 +8,6 @@ LIME=/p/home/jusers/boyle2/juwels/gm2dwf/boyle/
     --disable-gparity \
     --disable-fermion-reps \
     --with-lime=$LIME \
-    --enable-accelerator-cshift \
     --disable-unified \
     CXX=nvcc \
     LDFLAGS="-cudart shared " \

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

systems/Frontier/sourceme-rocm6.4.sh

@@ -0,0 +1,14 @@
+
+echo spack
+. /autofs/nccs-svm1_home1/paboyle/spack/share/spack/setup-env.sh
+
+export CLIME=`spack find --paths c-lime | grep ^c-lime | awk '{print $2}' `
+export MPFR=`spack find --paths mpfr    | grep ^mpfr  | awk '{print $2}' `
+export OPENSSL=`spack find --paths openssl | grep openssl | awk  '{print $2}' `
+export GMP=`spack find --paths gmp      | grep ^gmp | awk '{print $2}' `
+
+module load cce/20.0.0
+module load cpe/25.09
+module load rocm/6.4.0
+export LD_LIBRARY_PATH=$CRAY_LD_LIBRARY_PATH:$LD_LIBRARY_PATH
+export LD_LIBRARY_PATH=/opt/rocm-6.4.0/lib/llvm/lib/:$LD_LIBRARY_PATH

systems/Frontier/sourceme-rocm7.0.sh

@@ -0,0 +1,10 @@
+
+echo spack
+. /autofs/nccs-svm1_home1/paboyle/spack/share/spack/setup-env.sh
+
+
+module load cce/21.0.0
+module load cpe/26.03
+module load rocm/7.0.2
+export LD_LIBRARY_PATH=$CRAY_LD_LIBRARY_PATH:$LD_LIBRARY_PATH
+export LD_LIBRARY_PATH=/opt/rocm-7.0.2/lib/llvm/lib/:$LD_LIBRARY_PATH

systems/Frontier/sourceme-rocm7.2.sh

@@ -0,0 +1,14 @@
+
+echo spack
+. /autofs/nccs-svm1_home1/paboyle/spack/share/spack/setup-env.sh
+
+export CLIME=`spack find --paths c-lime | grep ^c-lime | awk '{print $2}' `
+export MPFR=`spack find --paths mpfr    | grep ^mpfr  | awk '{print $2}' `
+export OPENSSL=`spack find --paths openssl | grep openssl | awk  '{print $2}' `
+export GMP=`spack find --paths gmp      | grep ^gmp | awk '{print $2}' `
+
+module load cce/21.0.0
+module load cpe/26.03
+module load rocm/7.2.0
+export LD_LIBRARY_PATH=$CRAY_LD_LIBRARY_PATH:$LD_LIBRARY_PATH
+export LD_LIBRARY_PATH=/opt/rocm-7.2.0/lib/llvm/lib/:$LD_LIBRARY_PATH

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
+module load cce/21.0.0
-export LD_LIBRARY_PATH=/opt/cray/libfabric/1.20.1/lib64/:$LD_LIBRARY_PATH
+module load cpe/26.03
-export LD_LIBRARY_PATH=/opt/gcc/mpfr/3.1.4/lib:$LD_LIBRARY_PATH
+module load rocm/7.0.2
-export LD_LIBRARY_PATH=`pwd`/:$LD_LIBRARY_PATH
+export LD_LIBRARY_PATH=$CRAY_LD_LIBRARY_PATH:$LD_LIBRARY_PATH
-export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$HOME/LD_PATH/
+export LD_LIBRARY_PATH=/opt/rocm-7.0.2/lib/llvm/lib/:$LD_LIBRARY_PATH
-
-#echo spack load c-lime
-#spack load c-lime
-#module load emacs 
-##module load PrgEnv-gnu
-##module load cray-mpich
-##module load cray-fftw
-##module load craype-accel-amd-gfx90a
-##export LD_LIBRARY_PATH=/opt/gcc/mpfr/3.1.4/lib:$LD_LIBRARY_PATH
-#Hack for lib
-##export LD_LIBRARY_PATH=`pwd`/:$LD_LIBRARY_PATH

systems/Perlmutter/config-command

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

systems/Perlmutter/sourceme.sh

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

systems/mac-arm/config-command-mpi

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

systems/mac-arm/sourceme.sh

@@ -0,0 +1,11 @@
+source /Users/peterboyle/QCD//Spack/spack//share/spack/setup-env.sh
+
+export FFTW=`spack find --paths fftw    | grep ^fftw   | awk '{print $2}' `
+#export HDF5=`spack find --paths hdf5+cxx   | grep ^hdf5   | awk '{print $2}' `
+export CLIME=`spack find --paths c-lime | grep ^c-lime | awk '{print $2}' `
+export MPFR=`spack find --paths mpfr    | grep ^mpfr  | awk '{print $2}' `
+export OPENSSL=`spack find --paths openssl | grep openssl | awk  '{print $2}' `
+export GMP=`spack find --paths gmp      | grep ^gmp | awk '{print $2}' `
+
+export LD_LIBRARY_PATH=$MPFR/lib:$LD_LIBRARY_PATH
+export LD_LIBRARY_PATH=$GMP/lib:$LD_LIBRARY_PATH

tests/Test_fft_memory.cc

@@ -0,0 +1,261 @@
+/*************************************************************************************
+
+ Test_fft_memory.cc
+
+ Memory growth test for PlannedFFT on a spin-colour matrix (propagator) field.
+
+ The test creates a single PlannedFFT object (which allocates FFTW plans once),
+ then repeatedly applies FFT_all_dim to the same propagator 400 times.
+
+ If PlannedFFT is working correctly the RSS should remain flat after the first
+ iteration — no new plans, no new deviceVector allocations beyond the per-call
+ pencil buffer which is freed at the end of each FFT_dim_execute call.
+
+ Build exactly like any other Grid test, e.g.:
+   make Test_fft_memory
+ or compile manually:
+   $(CXX) $(CXXFLAGS) Test_fft_memory.cc -o Test_fft_memory $(LDFLAGS)
+
+*************************************************************************************/
+
+#include <Grid/Grid.h>
+using namespace Grid;
+
+// --------------------------------------------------------------------------
+// Helper: read RSS (resident set size) in kB from /proc/self/status.
+// Returns 0 on platforms where /proc is unavailable.
+// --------------------------------------------------------------------------
+static long getCPURSSKb()
+{
+  long rss = 0;
+  FILE *fp = fopen("/proc/self/status", "r");
+  if (!fp) return -1;
+  char line[256];
+  while (fgets(line, sizeof(line), fp)) {
+    if (strncmp(line, "VmRSS:", 6) == 0) {
+      sscanf(line + 6, "%ld", &rss);
+      break;
+    }
+  }
+  fclose(fp);
+  return rss;
+}
+
+static long getGPUUsedMb()
+{
+#if defined(GRID_CUDA)
+  size_t free_bytes  = 0;
+  size_t total_bytes = 0;
+  cudaError_t err = cudaMemGetInfo(&free_bytes, &total_bytes);
+  if (err != cudaSuccess) return -1;
+  return (long)((total_bytes - free_bytes) / (1024 * 1024));
+ 
+#elif defined(GRID_HIP)
+  size_t free_bytes  = 0;
+  size_t total_bytes = 0;
+  hipError_t err = hipMemGetInfo(&free_bytes, &total_bytes);
+  if (err != hipSuccess) return -1;
+  return (long)((total_bytes - free_bytes) / (1024 * 1024));
+ 
+#else
+  return -1;   // CPU-only build: no GPU to query
+#endif
+}
+
+// ============================================================
+//  Convenience struct — one snapshot of both sides
+// ============================================================
+struct MemSnapshot {
+  long cpu_rss_kb;   // host RSS in kB  (-1 if unavailable)
+  long gpu_used_mb;  // device used in MB (-1 if no GPU)
+};
+ 
+static MemSnapshot takeSnapshot()
+{
+  MemSnapshot s;
+  s.cpu_rss_kb  = getCPURSSKb();
+  s.gpu_used_mb = getGPUUsedMb();
+  return s;
+}
+ 
+// ============================================================
+//  Pretty-print one row of the monitoring table
+// ============================================================
+static void printRow(int iter,
+                     const MemSnapshot &now,
+                     const MemSnapshot &prev)
+{
+  long cpu_delta = (now.cpu_rss_kb  >= 0 && prev.cpu_rss_kb  >= 0)
+                 ? now.cpu_rss_kb  - prev.cpu_rss_kb  : 0;
+  long gpu_delta = (now.gpu_used_mb >= 0 && prev.gpu_used_mb >= 0)
+                 ? now.gpu_used_mb - prev.gpu_used_mb : 0;
+ 
+  // Sign prefix so deltas are unambiguous
+  auto sign = [](long v) -> const char* { return v >= 0 ? "+" : ""; };
+ 
+  std::cout << GridLogMessage
+            << std::setw(6)  << iter
+            << "  CPU: " << std::setw(10) << now.cpu_rss_kb  << " kB"
+            << "  (" << sign(cpu_delta) << std::setw(7) << cpu_delta << " kB)"
+            << "  GPU: " << std::setw(7)  << now.gpu_used_mb << " MB"
+            << "  (" << sign(gpu_delta) << std::setw(5) << gpu_delta << " MB)"
+            << "\n";
+}
+
+// ============================================================
+
+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;
+
+  // ------------------------------------------------------------------
+  // Grid setup — use whatever lattice/mpi/simd was passed on the CLI,
+  // e.g.  --grid 8.8.8.8 --mpi 1.1.1.1
+  // ------------------------------------------------------------------
+  Coordinate latt_size   = GridDefaultLatt();
+  Coordinate simd_layout = GridDefaultSimd(Nd, vComplexD::Nsimd());
+  Coordinate mpi_layout  = GridDefaultMpi();
+
+  GridCartesian GRID(latt_size, simd_layout, mpi_layout);
+
+  int vol = 1;
+  for (int d = 0; d < (int)latt_size.size(); d++) vol *= latt_size[d];
+
+  std::cout << GridLogMessage << "Lattice : ";
+  for (int d = 0; d < Nd; d++) std::cout << latt_size[d] << " ";
+  std::cout << std::endl;
+
+  // ------------------------------------------------------------------
+  // Propagator field: SpinColourMatrix = 12x12 complex, i.e.
+  // LatticePropagatorD (= Lattice<iSpinColourMatrix<vComplexD>>).
+  // This is the standard QCD quark propagator type.
+  // ------------------------------------------------------------------
+  LatticePropagatorD prop(&GRID);
+  
+  // ------------------------------------------------------------------
+  // Fill the propagator with a momentum-space plane wave,
+  // following the pattern from Test_fft.cc.
+  // We set each spin-colour component (a,b) to exp(i * sum_mu p_mu x_mu)
+  // with a fixed momentum p = (1,2,1,2).
+  // ------------------------------------------------------------------
+  Coordinate pvec({1, 2, 1, 2});
+
+  LatticeComplexD phase(&GRID);
+  LatticeComplexD coor(&GRID);
+  ComplexD ci(0.0, 1.0);
+
+  phase = Zero();
+  for (int mu = 0; mu < Nd; mu++) {
+    RealD TwoPiL = M_PI * 2.0 / latt_size[mu];
+    LatticeCoordinate(coor, mu);
+    phase = phase + (TwoPiL * pvec[mu]) * coor;
+  }
+  phase = exp(phase * ci);   // e^{i p.x}
+
+  // Broadcast the phase into every spin-colour matrix entry
+  prop = Zero();
+  prop = prop + phase;
+
+  std::cout << GridLogMessage
+            << "Propagator norm2 = " << norm2(prop) << std::endl;
+
+  // ------------------------------------------------------------------
+  // Baseline snapshot BEFORE PlannedFFT construction
+  // ------------------------------------------------------------------
+  MemSnapshot snap_before_plan = takeSnapshot();
+  std::cout << GridLogMessage
+            << "[mem] Before PlannedFFT construction"
+            << "  CPU: " << snap_before_plan.cpu_rss_kb  << " kB"
+            << "  GPU: " << snap_before_plan.gpu_used_mb << " MB"
+            << std::endl;
+ 
+  // ------------------------------------------------------------------
+  // Create the PlannedFFT — plans are allocated here ONCE for all
+  // dimensions and stored inside the object.
+  // ------------------------------------------------------------------
+  PlannedFFT<iSpinColourMatrix<vComplexD>> plannedFFT(&GRID);
+
+  // ------------------------------------------------------------------
+  // Snapshot AFTER plan construction — this is the true baseline
+  // for the loop, because cufftPlanMany itself grabs device memory.
+  // ------------------------------------------------------------------
+  MemSnapshot snap_after_plan = takeSnapshot();
+  std::cout << GridLogMessage
+            << "[mem] After  PlannedFFT construction"
+            << "  CPU: " << snap_after_plan.cpu_rss_kb  << " kB"
+            << "  GPU: " << snap_after_plan.gpu_used_mb << " MB"
+            << "  (plan overhead:"
+            << "  CPU +" << snap_after_plan.cpu_rss_kb - snap_before_plan.cpu_rss_kb << " kB"
+            << "  GPU +" << snap_after_plan.gpu_used_mb - snap_before_plan.gpu_used_mb << " MB)"
+            << std::endl;
+
+  MemoryManager::Print();
+  // ------------------------------------------------------------------
+  // 400-iteration loop.
+  // Each iteration computes the full 4d forward FFT of `prop`.
+  // We deliberately do NOT cache the result — we always start from
+  // the same `prop` so the FFT is recomputed identically each time.
+  // The point is to watch memory, not correctness.
+  // ------------------------------------------------------------------
+  const int Niter = 40;
+  const int Niter2 = 32;
+
+  // Print header for the memory table
+  std::cout << GridLogMessage
+              << "\n"
+              << std::setw(6)  << "iter"
+              << "  CPU: " << std::setw(10) << "RSS[kB]"
+              << "  (  delta  )"
+              << "  GPU: " << std::setw(7)  << "used[MB]"
+              << "  (delta)"
+              << "\n";
+ 
+  MemSnapshot snap_prev = snap_after_plan;
+
+  for (int i = 0; i < Niter; i++) {
+    std::vector<LatticePropagatorD> G;
+
+    for (int j = 0; j < Niter2; j++) {
+      LatticePropagatorD prop_fft(&GRID);
+
+      // Full 4d forward FFT using the pre-built plans
+      plannedFFT.FFT_all_dim(prop_fft, prop, FFT::forward);
+
+      G.push_back(prop_fft);
+    }
+ 
+    // cudaMemGetInfo reflects the state *after* any pooled frees have
+    // been committed, so this is accurate without an explicit sync —
+    // FFT_dim_execute already calls accelerator_barrier() internally.
+    MemSnapshot snap_now = takeSnapshot();
+    printRow(i, snap_now, snap_prev);
+    MemoryManager::Print();
+    snap_prev = snap_now;
+  }
+
+  // ------------------------------------------------------------------
+  // Summary
+  // ------------------------------------------------------------------
+  MemSnapshot snap_final = takeSnapshot();
+ 
+  long cpu_growth = snap_final.cpu_rss_kb  - snap_after_plan.cpu_rss_kb;
+  long gpu_growth = snap_final.gpu_used_mb - snap_after_plan.gpu_used_mb;
+ 
+  std::cout << GridLogMessage
+            << "\n==== Memory summary (baseline = after plan construction) ====\n"
+            << "  CPU RSS growth over " << Niter << " FFTs : "
+            << cpu_growth << " kB"
+            << (cpu_growth == 0 ? "  OK" : "  *** GROWING ***") << "\n"
+            << "  GPU used growth over " << Niter << " FFTs : "
+            << gpu_growth << " MB"
+            << (gpu_growth == 0 ? "  OK" : "  *** GROWING ***") << "\n"
+            << "  Note: first-call watermark from pool fill is expected and benign.\n"
+            << "        A leak shows as continuous growth beyond iter ~2-3.\n";
+ 
+  Grid_finalize();
+  return 0;
+}

tests/core/Test_fft.cc

@@ -113,7 +113,6 @@ 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,7 +95,6 @@ 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/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();
+}