1D порционный БПФ вещественных массивов

Question

Я слышал/читал, что мы можем использовать пакетный режим CUFFT, если у нас есть некоторые п БПФА для выполнения некоторых м векторов каждых. Поэтому, чтобы проверить это, я сделал образец программы и запустил ее. Данные, которые я использовал, были файлом с числами 1024 с плавающей запятой, так же, как и номера 1024, повторяющиеся 10 раз. В то время как я должен получить тот же результат для 1024 point 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

Как отметил Роберт Crovella, и, как сообщили в руководстве пользователя CUFFT - CUDA 6.5,

Пакетный размеры, кроме 1 для cufftPlan1d() было осуждается. Используйте команду 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); 

Имейте в виду, что ваш пример не учитывает тот факт, что 1D FFT из cufftReal массива длины 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 error handling.