CUDA cufft 2D пример

Question

В настоящее время я работаю над программой, которая должна реализовать 2D-FFT (для взаимной корреляции). Я сделал 1D FFT с CUDA, который дал мне правильные результаты, теперь я пытаюсь реализовать 2D-версию. С несколькими примерами и документацией онлайн мне трудно понять, что такое ошибка.

До сих пор я использовал только руководство cuFFT.

В любом случае, я создал два массива 5x5 и наполнил их 1-м. Я скопировал их на память GPU и сделал передовой FFT, умножил их и затем сделал ifft на результат. Это дает мне массив 5x5 со значениями 650. Я ожидаю получить сигнал постоянного тока со значением 25 только в одном слоте в массиве 5x5. Вместо этого я получаю 650 во всем массиве.

Кроме того, мне не разрешено распечатывать значение сигнала после его копирования в память графического процессора. Письмо

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 cufftExecR2C\n"); 
    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 cufftExecC2C\n"); 
    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 

Answer 1

Это дает мне массив 5x5 со значениями 650: Он читает 625, который 5 * 5 * 5 * 5. Алгоритм свертки, который вы используете, требует дополнительного деления на N * N. Действительно, в cufft коэффициент преобразования в прямом преобразовании отсутствует. Следовательно, ваша свертка не может быть простым умножением двух полей в частотной области. (некоторые называли бы это математиками DFT, а не врачами DFT).

Кроме того я не допускается, чтобы распечатать значение сигнала после того, как он был скопирован на память GPU: Это стандартное поведение CUDA. При распределении памяти на устройстве данные существуют в адресном пространстве памяти устройства и не могут быть доступны ЦП без дополнительных усилий. Поиск управляется памятью, или zerocopy для доступа к данным с обеих сторон PCI Express (об этом говорится во многих других сообщениях).

Answer 2

Есть несколько проблем здесь:

1. Вы запускаете слишком много потоков для размера входных массивов в ядре умножения, так что должно быть неудачей с недоступным ошибком памяти. Я удивлен, что вы не получаете какой-либо ошибки времени выполнения.
2. Ваше ожидаемое решение от fft/fft-dot-продукта - если это так, я считаю неверным. Правильным решением будет матрица 5x5 с 25 в каждой записи.
3. Как ясно описан в документации CUFFT, библиотека выполняет unnormalised: БПФ

   CUFFT выполняет уна-нормированный БПФ; то есть выполнение прямого БПФ по набору входных данных с последующим обратным БПФ на полученном наборе дает данные, которые равны входному сигналу, масштабируемому по количеству элементов. Масштабирование или преобразование по обратному размеру набора данных остается для пользователя, чтобы он выполнялся по мере возможности.

Так по моему расчету, правильное решение для вывода кода должен быть 5x5 матрица с 625 в каждой записи, которая будет нормированной матрицы 5х5 с 25 в каждой записи, то есть. ожидаемый результат. Я не понимаю, как проблема в (1) не приводит к разным результатам, поскольку ядро ​​умножения должно терпеть неудачу.

TLDR; нечего здесь видеть, двигаться дальше ...