GPU加速，利用CUFFT快速实现二维FFT运算

下面再次用一个简单的例子来介绍在c-mex函数中怎样调用CUFFT函数库。本例中，将随机生成一个二维实矩阵，进行二维FFT运算，并将其复数结果返回MATLAB。

1.步骤1

打开MATLAB命令窗口，创建一个新文件，将其保存为cufftDemo.cpp。

2.步骤2

在空文件cuftDemo.cpp中，输入以下代码：

与CUBLAS例子中相同，这是空子例行程序。

3.步骤3

在这步中，检查输入数据类型是单精度还是浮点型：

4.步骤4

这步中，我们在设备中为输入矩阵A创建存储空间，并将数据复制到GPU设备中：

5.步骤5

现在，在GPU中创建空间存储FFT输出的复数据

这里我们必须要格外注意，MATLAB数据是按列存储，这表示同一列中的数据在内存空间中是连续存放的。然而，CUFFT数据是按行存储的。CUFFT希望数据是按行连续存储，而不是像MATLAB那样按列存储。

表6.5是CUFFT API参考中给出的输入输出数据大小汇总。这里的2D和R2C（实数-复数）分别是我们关注的维度和类型。N₂是按列变换的数据，当在内存空间中连续移动时，N₂快速变化。所以在MATLAB数据类型中，N₁是行维度，N₂是列维度，MATLAB中的列维度和行维度在CUFFT中分别对应N₁和N₂。如果我们不注意维度和数据布局，得出的结果将发生转置。

表6.5 输入和输出数据大小（来自CUFFT AP参考

●MATLAB中的行维度M对应CUFFT中的N2。

●MATLAB中的列维度N对应CUFFT中的N1。

所以，输出的行数是实际行数的一半加1。看一下函数参考中给出的函数cufftPlan2d(…)的定义，这个函数将在下一步中调用。

cufftResult cufftPlan2d(cufftHandle *plan, int nx, int ny, cufftType type)是根据指定的信号大小和数据类型创建2D FFT配置。例如，输入信息如下：

对于输入参数nx和ny，我们分别传递numACols和numARows的值。(www.xing528.com)

6.步骤6

输入和输出数据矩阵准备就绪后，调用CUFFT函数。像CUBLAS例子中那样，首先创建它的句柄：

在cufftPlan2d(…)中，我们确定了输入矩阵的大小以及希望进行何种运算。一旦制定了这一计划，在之后的程序中如果需要的话可以再次利用，不需要时可以调用cufftDestroy(…)函数将其清除。实际的计算在cufftExecR2C(…)函数中进行，在输出数组deviceOut中可得到输出。

7.步骤7

在GPU内存中可以得到FFT计算数据。将这些计算数据返回到主机存储器中，之后就可以在MATLAB中使用它了。

与CUBLAS不同，我们用CUDA函数分配和释放内存，且在GPU和c-mex之间相互传送数据。但除了这三个CUDA函数外，不需要定义我们自己的内核函数，也不需要确定线程块和线程的大小。正如第4章中讨论的那样，MATLAB把复数存放在两个独立的矩阵中：一个存储实部，一个存储虚部。CUFFT将复数一个一个存储在一块内存空间中。所以在这步中，我们要为MATLAB创建一个复数矩阵，并把实部和虚部数据复制到相应的矩阵中。

8.步骤8

现在我们已经实现了c-mex代码。进入MATLAB命令窗口，然后编译我们的代码：

各选项含义如下：

●-lcudart：表明在调用CUDA运行时库。具体而言，我们调用了两个基本CUDA函数：cudaMalloc(…)和cudaFree(…)。

●-lcufft：表明在调用CUFFT库。具体而言，我们调用了cufftXXX函数调用。

●-Ldir：dir是CUDA和CUFFT库所在的目录。

●-Idir：dir是CUDA和CUFFT头文件所在的目录。

9.步骤9

如在64位Windows操作系统下进行编译，前面c-mex代码编译会产生c-mex文件cublasDemo.mexw64。在MATLAB命令窗口中调用我们的函数来进行以下乘法：

比较两个结果：一个用MATLAB fft2函数，另一个用CUFFT c-mex函数。注意到cufftDemo(A)返回结果比MATLAB fft2函数少一行，这一输出结果依据的是其API规范（见表6.5）。

本例中，N₁=4，N₂=4，这样最终输出结果为三行。

在这个例子中，缺少的这行可以用第二行的复共轭计算得到，如下所示：