stateful_fftw.cpp

struct plan_pkg {
  unsigned N1, N2, N3;
  ULONG N;
  ULONG Nhalf;
  complex_prec *C;
  real_prec *R;
  #ifdef SINGLE_PREC
  fftwf_plan plan;
  #endif
  #ifdef DOUBLE_PREC
  fftw_plan plan;
  #endif

  // R2C
  plan_pkg(unsigned int _N1, unsigned int _N2, unsigned int _N3, real_prec *in, complex_prec *out) {
    std::cout << "plan_pkg R2C" << std::endl;

    if (_N1 > INT_MAX || _N2 > INT_MAX || _N3 > INT_MAX) {
      std::string msg = "FFT dimensions must not be larger than INT_MAX!";
      throw std::logic_error(msg);
    }

    N1 = _N1, N2 = _N2, N3 = _N3;
    N = static_cast<ULONG>(N1*N2*N3);
    Nhalf = static_cast<ULONG>(N1 * N2 * (N3/2 + 1));
    R = in;
    C = out;

    #ifdef SINGLE_PREC
    plan = fftwf_plan_dft_r2c_3d(static_cast<int>(N1), static_cast<int>(N2), static_cast<int>(N3), in, out, FFTW_PATIENT);
    #endif
    #ifdef DOUBLE_PREC
    plan = fftw_plan_dft_r2c_3d(static_cast<int>(N1), static_cast<int>(N2), static_cast<int>(N3), in, out, FFTW_PATIENT);
    #endif
  }

  // C2R
  plan_pkg(unsigned int _N1, unsigned int _N2, unsigned int _N3, complex_prec *in, real_prec *out) {
    std::cout << "plan_pkg C2R" << std::endl;

    if (_N1 > INT_MAX || _N2 > INT_MAX || _N3 > INT_MAX) {
      std::string msg = "FFT dimensions must not be larger than INT_MAX!";
      throw std::logic_error(msg);
    }

    N1 = _N1, N2 = _N2, N3 = _N3;
    N = static_cast<ULONG>(N1*N2*N3);
    Nhalf = static_cast<ULONG>(N1 * N2 * (N3/2 + 1));
    C = in;
    R = out;

    #ifdef SINGLE_PREC
    plan = fftwf_plan_dft_c2r_3d(static_cast<int>(N1), static_cast<int>(N2), static_cast<int>(N3), in, out, FFTW_PATIENT);
    #endif
    #ifdef DOUBLE_PREC
    plan = fftw_plan_dft_c2r_3d(static_cast<int>(N1), static_cast<int>(N2), static_cast<int>(N3), in, out, FFTW_PATIENT);
    #endif
  }

  ~plan_pkg() {
    #ifdef SINGLE_PREC
    fftwf_destroy_plan(plan);
    #endif
    #ifdef DOUBLE_PREC
    fftw_destroy_plan(plan);
    #endif
  }
};


// Real FFT versions that take a plan as input. This greatly speeds up the FFT

// Take care: multidimensional C2R fftw transforms destroy the input array!
void fftC2Rplanned(complex_prec *in, real_prec *out, struct plan_pkg *plan) {
  if (in != plan->C) {
    copyArray(in, plan->C, plan->Nhalf);
  }

  #ifdef SINGLE_PREC
  fftwf_execute(plan->plan);
  #endif
  #ifdef DOUBLE_PREC
  fftw_execute(plan->plan);
  #endif

  real_prec fac = 1 / static_cast<real_prec>(plan->N);
  multiply_factor_array(fac, plan->R, out, plan->N);
}

void fftR2Cplanned(real_prec *in, complex_prec *out, struct plan_pkg *plan) {
  if (in != plan->R) {
    copyArray(in, plan->R, plan->N);
  }

  #ifdef SINGLE_PREC
  fftwf_execute(plan->plan);
  #endif
  #ifdef DOUBLE_PREC
  fftw_execute(plan->plan);
  #endif

  if (out != plan->C) {
    copyArray(plan->C, out, plan->Nhalf);
  }
}