基于cuFFT的傅里叶变换：复-复变换是否更高效？

我正在用CUDA/C++编写傅里叶空间中偏微分方程（PDE）的时间积分代码，需要演化一个实值数组。我实现了两种逻辑完全一致的方案：

方案一：全复数数组+复-复变换 #

所有数组均定义为复数类型，执行复-复傅里叶变换。为保证实值数组始终保持实值，每个时间步会将其虚部置0，消除累积的微小数值误差。

float *real;
cufftComplex *fourier;
...
for (int t = 0; t < totalTime; t++)
{
   cufftExecC2R(plan1, fourier, real);
   // 实空间计算操作
   cufftExecR2C(plan2, real, fourier);
   integrate<<<blocks, threads>>>(fourier, other parameters);
}

方案二：实数组+实-复/复-实变换 #

用浮点数组存储实值数据，复数数组存储傅里叶变换结果，执行复-实和实-复变换。理论上该方案内存效率更高：实数组替代复数组，傅里叶空间的复数组仅需存储N/2+1个元素。

cufftComplex *real;
cufftComplex *fourier;
...
for (int t = 0; t < totalTime; t++)
{
   cufftExecC2C(plan, fourier, real, CUFFT_INVERSE);
   // 实空间计算操作
   cufftExecC2C(plan, real, fourier, CUFFT_FORWARD); // 仅需一个变换计划
   integrate<<<blocks, threads>>>(fourier, other parameters);
}

但实际测试中，方案一在我的硬件上始终比方案二快约20%。这个结果符合预期吗？复-复变换在计算效率上真的优于实-复/复-实变换？

最小复现代码 #

以下两段代码可直接编译。在我的硬件上，方案二平均耗时1.1-1.2秒，方案一平均耗时约0.9秒。由于实际程序运行时间极长，这一微小差异会被显著放大。

helper_cuda.h和helper_functions.h并非编译必需文件，可在CUDA Toolkit示例代码的Common目录中找到。

仅使用复-复变换的代码 #

#include <iostream>
#include <fstream>
#include <math.h>
#include <random>
#include <fftw3.h>
#include <complex.h>


#include <cuda_runtime.h>
#include <cufft.h>
#include <cufftXt.h>
#include <helper_cuda.h>
#include <helper_functions.h>

using namespace std;

void printProgressbar(int t, int Nsteps, int length)
{
    int progress = t*100/Nsteps;
    int sprogress = t*length/Nsteps;
    int left = length - sprogress - 1;
    cout << "Progress [";
    for (int i = 0; i < sprogress; i++)
        cout << "=";
    cout << ">";
    for (int i = 0; i < left; i++)
        cout << " ";
    cout << "] " << progress << "%   \r";
    cout.flush();
}

int main()
{
    int N = 128;
    int Nsteps = 50000;
    int numberofsnaps = 5;
    int plot_steps = Nsteps/numberofsnaps;

    cufftComplex *phi_host = new cufftComplex[N*N];
    random_device rd;
    mt19937 gen(rd());
    uniform_real_distribution<> dis(-0.01,0.01);
    for (int i = 0; i < N*N; i++)
    {
        phi_host[i].x = (float)dis(gen);
        phi_host[i].y = 0.0f;
    }

    // Device arrays
    cufftComplex *phi_real;
    cufftComplex *phi_comp;

    cudaMalloc(reinterpret_cast<void **>(&phi_real), N * N * sizeof(cufftComplex));
    cudaMalloc(reinterpret_cast<void **>(&phi_comp), N * N * sizeof(cufftComplex));

    cufftHandle plan_c2c;
    cufftPlan2d(&plan_c2c, N, N, CUFFT_C2C);

    cudaMemcpy(phi_real, phi_host, N * N * sizeof(cufftComplex), cudaMemcpyHostToDevice);
    cufftExecC2C(plan_c2c, phi_real, phi_comp, CUFFT_FORWARD);

    for (int t = 0; t < Nsteps; t++)
    {
        // Output data
        if (t % plot_steps == 0)
        {
            printProgressbar(t,Nsteps,25);
            cudaMemcpy(phi_host, phi_real, N * N * sizeof(cufftComplex), cudaMemcpyDeviceToHost);
        }
        cufftExecC2C(plan_c2c, phi_comp, phi_real, CUFFT_INVERSE);
        cufftExecC2C(plan_c2c, phi_real, phi_comp, CUFFT_FORWARD);
    }
    cout << "Progress [========================>] 100%" << endl;
    return 0;
}

使用实-复和复-实变换的代码 #

#include <iostream>
#include <fstream>
#include <math.h>
#include <random>
#include <fftw3.h>
#include <complex.h>


#include <cuda_runtime.h>
#include <cufft.h>
#include <cufftXt.h>
#include <helper_cuda.h>
#include <helper_functions.h>
using namespace std;

void printProgressbar(int t, int Nsteps, int length)
{
    int progress = t*100/Nsteps;
    int sprogress = t*length/Nsteps;
    int left = length - sprogress - 1;
    cout << "Progress [";
    for (int i = 0; i < sprogress; i++)
        cout << "=";
    cout << ">";
    for (int i = 0; i < left; i++)
        cout << " ";
    cout << "] " << progress << "%   \r";
    cout.flush();
}

int main()
{
    int N = 128;
    int Np = N/2 + 1;
    int Nsteps = 50000;
    int numberofsnaps = 5;
    int plot_steps = Nsteps/numberofsnaps;

    float *phi_host = new float[N*N];
    random_device rd;
    mt19937 gen(rd());
    uniform_real_distribution<> dis(-0.01,0.01);
    for (int i = 0; i < N*N; i++)
    {
        phi_host[i] = (float)dis(gen);
    }

    // Device arrays
    float *phi_real;
    cufftComplex *phi_comp;

    cudaMalloc(reinterpret_cast<void **>(&phi_real), N * N * sizeof(float));
    cudaMalloc(reinterpret_cast<void **>(&phi_comp), N * Np * sizeof(cufftComplex));

    cufftHandle plan_r2c, plan_c2r;
    cufftPlan2d(&plan_r2c, N, N, CUFFT_R2C);
    cufftPlan2d(&plan_c2r, N, N, CUFFT_C2R);

    cudaMemcpy(phi_real, phi_host, N * N * sizeof(float), cudaMemcpyHostToDevice);
    cufftExecR2C(plan_r2c, phi_real, phi_comp);

    for (int t = 0; t < Nsteps; t++)
    {
        // Output data
        if (t % plot_steps == 0)
        {
            printProgressbar(t,Nsteps,25);
            cudaMemcpy(phi_host, phi_real, N * N * sizeof(float), cudaMemcpyDeviceToHost);
        }
        cufftExecC2R(plan_c2r, phi_comp, phi_real);
        cufftExecR2C(plan_r2c, phi_real, phi_comp);
    }
    cout << "Progress [========================>] 100%" << endl;
    return 0;
}

内容的提问来源于stack exchange，提问作者MyUserIsThis