CUDA 袖带 2D 示例答案

【问题标题】：CUDA cufft 2D exampleCUDA 袖带 2D 示例
【发布时间】：2016-08-21 17:20:06
【问题描述】：

我目前正在开发一个必须实现 2D-FFT 的程序（用于互相关）。我用 CUDA 做了一个 1D FFT，它给了我正确的结果，我现在正在尝试实现一个 2D 版本。由于在线示例和文档很少，我很难找出错误所在。

到目前为止，我一直只使用 cuFFT 手册。

无论如何，我创建了两个 5x5 数组并用 1 填充它们。我已将它们复制到 GPU 内存并进行了前向 FFT，将它们相乘，然后对结果进行了 ifft。这给了我一个值为 650 的 5x5 阵列。我希望在 5x5 阵列的一个插槽中获得值为 25 的 DC 信号。相反，我在整个数组中得到 650。

此外，在将信号复制到 GPU 内存后，我不允许打印信号的值。写作

cout << d_signal[1].x << endl;

给我一个访问冲突。我在其他 cuda 程序中做过同样的事情，这不是问题。是否与复杂变量的工作方式有关，还是人为错误？

如果有人对出了什么问题有任何指示，我将不胜感激。这是代码

   #include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <helper_functions.h>
#include <helper_cuda.h>

#include <ctime>
#include <time.h>
#include <stdio.h>
#include <iostream>
#include <math.h>
#include <cufft.h>
#include <fstream>

using namespace std;
typedef float2 Complex;





__global__ void ComplexMUL(Complex *a, Complex *b)
{
    int i = threadIdx.x;
    a[i].x = a[i].x * b[i].x - a[i].y*b[i].y;
    a[i].y = a[i].x * b[i].y + a[i].y*b[i].x;
}


int main()
{


    int N = 5;
    int SIZE = N*N;


    Complex *fg = new Complex[SIZE];
    for (int i = 0; i < SIZE; i++){
        fg[i].x = 1; 
        fg[i].y = 0;
    }
    Complex *fig = new Complex[SIZE];
    for (int i = 0; i < SIZE; i++){
        fig[i].x = 1; // 
        fig[i].y = 0;
    }
    for (int i = 0; i < 24; i=i+5)
    {
        cout << fg[i].x << " " << fg[i + 1].x << " " << fg[i + 2].x << " " << fg[i + 3].x << " " << fg[i + 4].x << endl;
    }
    cout << "----------------" << endl;
    for (int i = 0; i < 24; i = i + 5)
    {
        cout << fig[i].x << " " << fig[i + 1].x << " " << fig[i + 2].x << " " << fig[i + 3].x << " " << fig[i + 4].x << endl;
    }
    cout << "----------------" << endl;

    int mem_size = sizeof(Complex)* SIZE;


    cufftComplex *d_signal;
    checkCudaErrors(cudaMalloc((void **) &d_signal, mem_size)); 
    checkCudaErrors(cudaMemcpy(d_signal, fg, mem_size, cudaMemcpyHostToDevice));

    cufftComplex *d_filter_kernel;
    checkCudaErrors(cudaMalloc((void **)&d_filter_kernel, mem_size));
    checkCudaErrors(cudaMemcpy(d_filter_kernel, fig, mem_size, cudaMemcpyHostToDevice));

    // cout << d_signal[1].x << endl;
    // CUFFT plan
    cufftHandle plan;
    cufftPlan2d(&plan, N, N, CUFFT_C2C);

    // Transform signal and filter
    printf("Transforming signal cufftExecR2C\n");
    cufftExecC2C(plan, (cufftComplex *)d_signal, (cufftComplex *)d_signal, CUFFT_FORWARD);
    cufftExecC2C(plan, (cufftComplex *)d_filter_kernel, (cufftComplex *)d_filter_kernel, CUFFT_FORWARD);

    printf("Launching Complex multiplication<<< >>>\n");
    ComplexMUL <<< 32, 256 >> >(d_signal, d_filter_kernel);

    // Transform signal back
    printf("Transforming signal back cufftExecC2C\n");
    cufftExecC2C(plan, (cufftComplex *)d_signal, (cufftComplex *)d_signal, CUFFT_INVERSE);

    Complex *result = new Complex[SIZE];
    cudaMemcpy(result, d_signal, sizeof(Complex)*SIZE, cudaMemcpyDeviceToHost);

    for (int i = 0; i < SIZE; i=i+5)
    {
        cout << result[i].x << " " << result[i + 1].x << " " << result[i + 2].x << " " << result[i + 3].x << " " << result[i + 4].x << endl;
    }

    delete result, fg, fig;
    cufftDestroy(plan);
    //cufftDestroy(plan2);
    cudaFree(d_signal);
    cudaFree(d_filter_kernel);

}

以上代码给出以下终端输出：

1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
----------------
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
----------------
Transforming signal cufftExecR2C
Launching Complex multiplication<<< >>>
Transforming signal back cufftExecC2C

625 625 625 625 625
625 625 625 625 625
625 625 625 625 625
625 625 625 625 625
625 625 625 625 625

【问题讨论】：

您发布的代码不完整，无法编译。你能解决这个问题吗？如果不编译和运行代码，很难告诉你可能出了什么问题，我现在做不到-
当然，我有一些我不想包含的未注释部分。我已将其删除并将所有内容编辑到我的帖子中。

标签： c++ cuda cufft

【解决方案1】：

这给了我一个 5x5 数组，其值为 650：它读取 625，即 5555。您使用的卷积算法需要除以 NN。确实，在cufft 中，前向变换中没有归一化系数。因此，您的卷积不能是频域中两个字段的简单相乘。（有些人会称之为数学家 DFT 而不是物理学家 DFT）。

此外，在将信号复制到 GPU 内存后，我不允许打印信号的值：这是标准 CUDA 行为。在设备上分配内存时，数据存在于设备内存地址空间中，不经额外努力就无法被 CPU 访问。搜索 managed 内存或 zerocopy 以从 PCI Express 的两侧访问数据（这在许多其他帖子中都有讨论）。

【讨论】：

感谢您的回复弗洛伦蒂感谢它。这对我帮助很大！

【解决方案2】：

这里有几个问题：

对于乘法内核中输入数组的大小而言，您启动了太多线程，因此应该会因内存越界错误而失败。我很惊讶您没有收到任何类型的运行时错误。
我相信您对 fft/fft - 点积 - ifft 序列的预期解决方案是不正确的。正确的解决方案是一个 5x5 矩阵，每个条目有 25 个。
正如 cuFFT 文档中明确描述的那样，该库执行 非标准化 FFT：

cuFFT 执行未归一化的 FFT；也就是说，对输入数据集执行正向 FFT，然后对结果集执行反向 FFT，得到的数据等于输入，按元素数量缩放。通过数据集大小的倒数缩放任一变换，留给用户执行。

因此，根据我的估计，您的代码的正确输出解决方案应该是一个 5x5 矩阵，每个条目中有 625，这将被归一化为每个条目中有 25 个的 5x5 矩阵，即。预期的结果。我不明白 (1) 处的问题如何没有产生不同的结果，因为乘法内核应该失败了。

TLDR;这里没什么可看的，继续前进……

【讨论】：

访问 GPU 上的有效内存区域，即使未分配，也不一定会在 cuda 内存检查测试之外发出错误。内核不一定会因为小溢出而失败。不过，您的所有积分仍然有效。
@FlorentDUGUET：输入数组是 25 个双字。内核启动每个块使用 256 个线程。当我运行它时（是的，我确实运行了它），它在 cuda-memcheck 中产生了数百个无效的内存访问错误。
感谢您的回复 talonmies 我很感激。这对我帮助很大！
@talonmies 我相信你有数百个无效的内存访问错误（我说不一定会在 cuda mem check 之外发出错误）。然而，内核并没有失败。我不知道这是否与将 GPU 共享到允许 CUDA 的多个 VM 的当前状态有关，而不仅仅是 OpenGL。

【解决方案3】：

就像已经提到的其他事情一样：我认为您的复数乘法内核没有做正确的事情。您将覆盖第一行中的a[i].x，然后使用第二行中a[i].x 的新值来计算a[i].y。我认为您需要在覆盖之前先生成a[i].x 的备份，例如：

float aReal_bk = a[i].x;
a[i].x = a[i].x * b[i].x - a[i].y * b[i].y;
a[i].y = aReal_bk * b[i].y + a[i].y * b[i].x;

【讨论】：