真实数组的一维批处理 FFT

Question

我听说/读过，如果我们有一些n FFT 来执行一些m 个向量，我们可以使用 cuFFT 的批处理模式。所以为了测试它，我制作了一个示例程序并运行它。我使用的数据是一个文件，其中包含一些1024浮点数，因为相同的1024数字重复了10多次。虽然1024点 FFT应该得到相同的结果，但我没有得到。如果我在某处概念上有错误，请纠正我，下面是代码，如果你能纠正我犯的一些错误。

注意：我只使用 1D FFT。

这是代码片段：

#include <cuda.h>
#include <cufft.h>
#include <stdio.h>
#include <math.h>

#define NX 1024
#define DATASIZE 1024
#define BATCH 10

int main (int argc, char* argv[])
{
        cufftHandle plan;
        cufftComplex *deviceOutputData, *hostOutputData;
        cufftReal *hostInputData, *deviceInputData;
        int i,j;

        FILE *in; // *out, *fp;

        cudaMalloc ((void**)&deviceInputData, NX*BATCH*sizeof(cufftReal));
        hostInputData = (cufftReal*) malloc (NX*BATCH*sizeof(cufftReal));

        cudaMalloc ((void**)&deviceOutputData, NX*BATCH*sizeof(cufftComplex));
        hostOutputData = (cufftComplex*) malloc (NX*BATCH*sizeof(cufftComplex));

        in = fopen ("InFile.txt", "r");

        if (in==NULL)
        {       fprintf (stderr, "Input file has some issues. Please check."); exit(1);}

        float data;
        //Allocate data

 for (i=0; i<BATCH; i++){
                for (j=0; j<DATASIZE;j++)
                {
                        fscanf(in, "%f", &data);
                        hostInputData [j + i*DATASIZE] = data;
                }
        }
        fclose (in);
        cudaMemcpy (deviceInputData, hostInputData, DATASIZE*BATCH*sizeof(cufftReal), cudaMemcpyHostToDevice);
        cufftPlan1d (&plan, NX, CUFFT_R2C, BATCH);
        cufftExecR2C (plan,  deviceInputData, deviceOutputData);
        cudaThreadSynchronize();
        cudaMemcpy (hostOutputData, deviceOutputData, DATASIZE*BATCH*sizeof(cufftComplex), cudaMemcpyDeviceToHost);
        cufftDestroy (plan);
        cudaFree (deviceOutputData);
        cudaFree (deviceInputData);

        #define a hostOutputData[j+i*NX].x
        #define b hostOutputData[j+i*NX].y
        float result[NX];
        for (i=0; i<BATCH; i++){
                printf ("\n*New Batch*\n");
                for (j=0; j<=NX/2;j++){
                        result[j] = sqrt ((a*a)+(b*b));
                        printf ("%f\n", result[j]);
                }

                for (j=1; j<NX/2; j++){
                        result[j+(NX/2)] = result [(NX/2)-j];
                        printf ("%f\n", result[j+(NX/2)]);
                }
        }

Answer 1

正如 Robert Crovella 所提到的，以及在 cuFFT 用户指南 - CUDA 6.5 中所报告的，

cufftPlan1d()不推荐使用1 以外的批次大小。使用cufftPlanMany()多个批处理执行。

下面，我报告了一个完整的示例，它更正了您的代码并使用cufftPlanMany()了cufftPlan1d(). 正如你将看到的，

int rank = 1;                           // --- 1D FFTs
int n[] = { DATASIZE };                 // --- Size of the Fourier transform
int istride = 1, ostride = 1;           // --- Distance between two successive input/output elements
int idist = DATASIZE, odist = (DATASIZE / 2 + 1); // --- Distance between batches
int inembed[] = { 0 };                  // --- Input size with pitch (ignored for 1D transforms)
int onembed[] = { 0 };                  // --- Output size with pitch (ignored for 1D transforms)
int batch = BATCH;                      // --- Number of batched executions
cufftPlanMany(&handle, rank, n, 
              inembed, istride, idist,
              onembed, ostride, odist, CUFFT_R2C, batch);

完全等同于“老式”

cufftPlan1d(&handle, DATASIZE, CUFFT_R2C, BATCH);

但是要注意，你的例子不占的事实的一维FFTcufftReal长度的数组DATASIZE是一个cufftComplex数组DATASIZE/2 + 1元素。

这是完整的示例：

#include <cuda.h>
#include <cufft.h>
#include <stdio.h>
#include <math.h>

#define DATASIZE 8
#define BATCH 2

/********************/
/* CUDA ERROR CHECK */
/********************/
#define gpuErrchk(ans) { gpuAssert((ans), __FILE__, __LINE__); }
inline void gpuAssert(cudaError_t code, const char *file, int line, bool abort=true)
{
   if (code != cudaSuccess) 
   {
      fprintf(stderr,"GPUassert: %s %s %d\n", cudaGetErrorString(code), file, line);
      if (abort) exit(code);
   }
}

/********/
/* MAIN */
/********/
int main ()
{
    // --- Host side input data allocation and initialization
    cufftReal *hostInputData = (cufftReal*)malloc(DATASIZE*BATCH*sizeof(cufftReal));
    for (int i=0; i<BATCH; i++)
        for (int j=0; j<DATASIZE; j++) hostInputData[i*DATASIZE + j] = (cufftReal)(i + 1);

    // --- Device side input data allocation and initialization
    cufftReal *deviceInputData; gpuErrchk(cudaMalloc((void**)&deviceInputData, DATASIZE * BATCH * sizeof(cufftReal)));
    cudaMemcpy(deviceInputData, hostInputData, DATASIZE * BATCH * sizeof(cufftReal), cudaMemcpyHostToDevice);

    // --- Host side output data allocation
    cufftComplex *hostOutputData = (cufftComplex*)malloc((DATASIZE / 2 + 1) * BATCH * sizeof(cufftComplex));

    // --- Device side output data allocation
    cufftComplex *deviceOutputData; gpuErrchk(cudaMalloc((void**)&deviceOutputData, (DATASIZE / 2 + 1) * BATCH * sizeof(cufftComplex)));

    // --- Batched 1D FFTs
    cufftHandle handle;
    int rank = 1;                           // --- 1D FFTs
    int n[] = { DATASIZE };                 // --- Size of the Fourier transform
    int istride = 1, ostride = 1;           // --- Distance between two successive input/output elements
    int idist = DATASIZE, odist = (DATASIZE / 2 + 1); // --- Distance between batches
    int inembed[] = { 0 };                  // --- Input size with pitch (ignored for 1D transforms)
    int onembed[] = { 0 };                  // --- Output size with pitch (ignored for 1D transforms)
    int batch = BATCH;                      // --- Number of batched executions
    cufftPlanMany(&handle, rank, n, 
                  inembed, istride, idist,
                  onembed, ostride, odist, CUFFT_R2C, batch);

    //cufftPlan1d(&handle, DATASIZE, CUFFT_R2C, BATCH);
    cufftExecR2C(handle,  deviceInputData, deviceOutputData);

    // --- Device->Host copy of the results
    gpuErrchk(cudaMemcpy(hostOutputData, deviceOutputData, (DATASIZE / 2 + 1) * BATCH * sizeof(cufftComplex), cudaMemcpyDeviceToHost));

    for (int i=0; i<BATCH; i++)
        for (int j=0; j<(DATASIZE / 2 + 1); j++)
            printf("%i %i %f %f\n", i, j, hostOutputData[i*(DATASIZE / 2 + 1) + j].x, hostOutputData[i*(DATASIZE / 2 + 1) + j].y);

    cufftDestroy(handle);
    gpuErrchk(cudaFree(deviceOutputData));
    gpuErrchk(cudaFree(deviceInputData));

}

请根据CUFFT 错误处理添加您自己的 cuFFT 错误检查。