Test_planned_fft: fix PlannedFFT template parameter to use ::vector_object

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

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

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