1D БПФ столбцов и строк 3D-матрицы в CUDA

Question

Я пытаюсь вычислить пакетные однофазные БПФ с использованием cufftPlanMany. Набор данных поступает из трехмерного поля, хранящегося в 1D-массиве, где я хочу вычислить 1D БПФ в направлениях x и y. Данные сохраняются, как показано на рисунке ниже; непрерывный в x, затем y, затем z.

Выполнение пакетных БПФ в x -направление (я считаю) несложное; с вводом stride=1, distance=nx и batch=ny * nz, он вычисляет БПФ над элементами {0,1,2,3}, {4,5,6,7}, ..., {28,29,30,31}. Тем не менее, я не могу придумать способ добиться того же для БПФ в направлении y. Партия для каждой плоскости xy снова является простой (ввод stride=nx, dist=1, batch=nx приводит к БПФ более {0,4,8,12}, {1,5,9,13} и т. Д.). Но с batch=nx * nz, перейдя с {3,7,11,15} до {16,20,24,28}, больше чем 1. Можно ли это сделать с помощью cufftPlanMany?

Answer 1

Я думаю, что короткий ответ на ваш вопрос (возможность использования одного cufftPlanMany для выполнения 1D БПФ из столбцов матрицы 3D) НЕТ.

Действительно, преобразования, выполненные в соответствии с cufftPlanMany, что вы называете, как

cufftPlanMany(&handle, rank, n, 
       inembed, istride, idist, 
       onembed, ostride, odist, CUFFT_C2C, batch); 

должны подчиняться Advanced Data Layout. В частности, 1D БПФ разработаны в соответствии со следующей компоновке

input[b * idist + x * istride] 

, где b обращается к b -го сигнала и istride является расстояние между двумя последовательными элементами в одном сигнале. Если 3D-матрица имеет размеры M * N * Q, и если вы хотите выполнить одномерные преобразования вдоль столбцов, то расстояние между двумя последовательными элементами будет M, а расстояние между двумя последовательными сигналами будет 1. Кроме того, количество партийных исполнений должно быть установлено равным M. С этими параметрами вы можете охватить только один фрагмент 3D-матрицы. Действительно, если вы попытаетесь увеличить M, тогда cuFFT начнет пытаться вычислить новые столбчатые БПФ, начиная со второй строки. Единственным решением этой проблемы является итеративный вызов cufftExecC2C для покрытия всех фрагментов Q.

Для записи следующий код предоставляет полностью обработанный пример того, как выполнять 1D БПФ столбцов 3D-матрицы.

#include <thrust/device_vector.h> 
#include <cufft.h> 

/********************/ 
/* 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); 
    } 
} 

int main() { 

    const int M = 3; 
    const int N = 4; 
    const int Q = 2; 

    thrust::host_vector<float2> h_matrix(M * N * Q); 

    for (int k=0; k<Q; k++) 
     for (int j=0; j<N; j++) 
      for (int i=0; i<M; i++) { 
       float2 temp; 
       temp.x = (float)(j + k * M); 
       //temp.x = 1.f; 
       temp.y = 0.f; 
       h_matrix[k*M*N+j*M+i] = temp; 
       printf("%i %i %i %f %f\n", i, j, k, temp.x, temp.y); 
      } 
    printf("\n"); 

    thrust::device_vector<float2> d_matrix(h_matrix); 

    thrust::device_vector<float2> d_matrix_out(M * N * Q); 

    // --- Advanced data layout 
    //  input[b * idist + x * istride] 
    //  output[b * odist + x * ostride] 
    //  b = signal number 
    //  x = element of the b-th signal 

    cufftHandle handle; 
    int rank = 1;       // --- 1D FFTs 
    int n[] = { N };      // --- Size of the Fourier transform 
    int istride = M, ostride = M;   // --- Distance between two successive input/output elements 
    int idist = 1, odist = 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 = M;       // --- Number of batched executions 
    cufftPlanMany(&handle, rank, n, 
        inembed, istride, idist, 
        onembed, ostride, odist, CUFFT_C2C, batch); 

    for (int k=0; k<Q; k++) 
     cufftExecC2C(handle, (cufftComplex*)(thrust::raw_pointer_cast(d_matrix.data()) + k * M * N), (cufftComplex*)(thrust::raw_pointer_cast(d_matrix_out.data()) + k * M * N), CUFFT_FORWARD); 
    cufftDestroy(handle); 

    for (int k=0; k<Q; k++) 
     for (int j=0; j<N; j++) 
      for (int i=0; i<M; i++) { 
       float2 temp = d_matrix_out[k*M*N+j*M+i]; 
       printf("%i %i %i %f %f\n", i, j, k, temp.x, temp.y); 
      } 

} 

Ситуация различна для случая, когда вы хотите выполнить 1D-преобразования строк. В этом случае расстояние между двумя последовательными элементами равно 1, тогда как расстояние между двумя последовательными сигналами равно M. Это позволяет вам установить число преобразований N * Q, а затем вызвать cufftExecC2C только один раз. Для записи приведенный ниже код дает полный пример 1D преобразований строк трехмерной матрицы.

#include <thrust/device_vector.h> 
#include <cufft.h> 

/********************/ 
/* 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); 
    } 
} 

int main() { 

    const int M = 3; 
    const int N = 4; 
    const int Q = 2; 

    thrust::host_vector<float2> h_matrix(M * N * Q); 

    for (int k=0; k<Q; k++) 
     for (int j=0; j<N; j++) 
      for (int i=0; i<M; i++) { 
       float2 temp; 
       temp.x = (float)(j + k * M); 
       //temp.x = 1.f; 
       temp.y = 0.f; 
       h_matrix[k*M*N+j*M+i] = temp; 
       printf("%i %i %i %f %f\n", i, j, k, temp.x, temp.y); 
      } 
    printf("\n"); 

    thrust::device_vector<float2> d_matrix(h_matrix); 

    thrust::device_vector<float2> d_matrix_out(M * N * Q); 

    // --- Advanced data layout 
    //  input[b * idist + x * istride] 
    //  output[b * odist + x * ostride] 
    //  b = signal number 
    //  x = element of the b-th signal 

    cufftHandle handle; 
    int rank = 1;       // --- 1D FFTs 
    int n[] = { M };      // --- Size of the Fourier transform 
    int istride = 1, ostride = 1;   // --- Distance between two successive input/output elements 
    int idist = M, odist = M;    // --- 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 = N * Q;      // --- Number of batched executions 
    cufftPlanMany(&handle, rank, n, 
        inembed, istride, idist, 
        onembed, ostride, odist, CUFFT_C2C, batch); 

    cufftExecC2C(handle, (cufftComplex*)(thrust::raw_pointer_cast(d_matrix.data())), (cufftComplex*)(thrust::raw_pointer_cast(d_matrix_out.data())), CUFFT_FORWARD); 
    cufftDestroy(handle); 

    for (int k=0; k<Q; k++) 
     for (int j=0; j<N; j++) 
      for (int i=0; i<M; i++) { 
       float2 temp = d_matrix_out[k*M*N+j*M+i]; 
       printf("%i %i %i %f %f\n", i, j, k, temp.x, temp.y); 
      } 

} 

Answer 2

Я предполагаю, что idist = nx * nz также мог бы перепрыгнуть целую плоскость, а пакет = nz затем покрыл бы одну плоскость yx. Решение должно быть принято в зависимости от того, больше ли nx или nz.