CuFFT Regression

example.cu

#include <cmath>

#include <algorithm>
#include <iostream>
#include <complex>

#include <cufft.h>

static constexpr int N1 = 64;
static constexpr int N2 = 512;

int main() {
  cufftComplex input_data[N2];
  for (int k = 0; k < N1; ++k) {
    input_data[k].x = 1.0f;
    input_data[k].y = 0.0f;
  }
  for (int k = N1; k < N2; ++k) {
    input_data[k].x = 0.0f;
    input_data[k].y = 0.0f;
  }

  cufftComplex* gpuData;
  cudaMalloc(&gpuData, N2*sizeof(*gpuData));
  cudaMemcpy(gpuData, input_data, N2*sizeof(*gpuData),
             cudaMemcpyHostToDevice);

  cufftHandle plan;
  cufftPlan1d(&plan, N2, CUFFT_C2C, 1);
  cufftExecC2C(plan, gpuData, gpuData, CUFFT_FORWARD);
  cufftDestroy(plan);

  cufftComplex output_data[N2];
  cudaMemcpy(output_data, gpuData, N2*sizeof(*gpuData),
             cudaMemcpyDeviceToHost);
  cudaFree(gpuData);

  double maxDiff = 0.0;
  for (int k = 0; k < N2; ++k) {
    double x = (M_PI*k)/N2;
    std::complex<double> z(0.0, (1 - N1)*x);
    z = std::exp(z);
    if (k > 0) z *= std::sin(N1*x) / std::sin(x);
    else z *= double(N1);

    std::complex<double> w(double(output_data[k].x),
                           double(output_data[k].y));
    maxDiff = std::max(std::fabs(z - w), maxDiff);
  }
  std::cout << "maxDiff=" << maxDiff << std::endl;
}