CUDA 袖口 2D 示例

CUDA cufft 2D example

本文关键字:示例 2D 袖口 CUDA      更新时间:2023-10-16

我目前正在开发一个必须实现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 cufftExecR2Cn");
    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 cufftExecC2Cn");
    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

这给了我一个值为 650 的 5x5 数组:它读取 625 ,即 5 555。您使用的卷积算法需要补充除以 N N。实际上,在袖口中,正向变换中没有归一化系数。因此,卷积不能是频域中两个场的简单乘法。(有些人会称它为数学家DFT,而不是物理学家DFT)。

此外,我不允许在将信号复制到 GPU 内存后打印出信号的值 :这是标准的 CUDA 行为。在设备上分配内存时,数据存在于设备内存地址空间中,CPU 无法在不增加工作的情况下访问这些数据。搜索托管内存或零拷贝以从PCI Express的两端访问数据(这在许多其他帖子中都有讨论)。

这里有几个问题:

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

    cuFFT 执行非规范化 FFT;也就是说,对输入数据集执行前向 FFT,然后在结果集上执行反向 FFT,生成的数据等于输入,按元素数缩放。按数据集大小的倒数缩放任一转换留给用户执行所需的操作。

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

TLDR;这里没什么可看的,继续前进...

就像已经提到的其他事情一样:我认为您的复杂乘法内核没有做正确的事情。您正在覆盖第一行中的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;