CUDA - Простое комплексное умножение чисел

Question

Я пытаюсь изучить CUDA, написав базовый код, который, надеюсь, поставит меня в лучшем положении для преобразования моего существующего кода на C++ в CUDA (для исследования).

Мне нужно сделать довольно сложную обработку чисел, поэтому я написал этот очень простой код, чтобы умножить массив комплексных чисел с реальным числом в ядре графического ядра.

#include <complex> 
#include <iostream> 
#include <cmath> 
#include "cuda.h" 
#include "math.h" 
#include "cuComplex.h" 

#define n 5 

using namespace std; 

#define gpuErrchk(ans) { gpuAssert((ans), __FILE__, __LINE__); } 
inline void gpuAssert(cudaError_t code, 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); 
    } 
} 

__global__ void func(double *s, cuDoubleComplex *j, cuDoubleComplex *calc) { 

    int tid = blockIdx.x; 

    calc[tid] = cuCmul(j[tid], make_cuDoubleComplex(*s, 0)); 

} 

int main(void) { 


    cuDoubleComplex calc[n+1], *dev_j, *dev_calc; 
    double *dev_s, s[n+1] = { 2.0, 2.0, 2.0, 2.0, 2.0 }; 
    //complex<double> j[n+1] 
    cuDoubleComplex j[n+1]; 

    for (int i = 1; i <= n; i++) { 
     j[i] = make_cuDoubleComplex(0, 5); 
     cout << "\nJ cout = " << cuCreal(j[i]) << ", " << cuCimag(j[i]); 
    } 

    // allocate the memory on the GPU 
    cudaMalloc((void**)&dev_s, (n+1) * sizeof(double)); 
    cudaMalloc((void**)&dev_j, (n+1) * sizeof(double)); 
    cudaMalloc((void**)&dev_calc, (n+1) * sizeof(double)); 

    cudaMemcpy(dev_s, s, (n+1) * sizeof(double), cudaMemcpyHostToDevice); 
    cudaMemcpy(dev_j, j, (n+1) * sizeof(double), cudaMemcpyHostToDevice); 

    func<<<n,1>>>(dev_s, dev_j, dev_calc); 
    //kernel<<<1,1>>>(a_d); 
    gpuErrchk(cudaPeekAtLastError()); 
    gpuErrchk(cudaMemcpy(calc, dev_calc, (n+1) * sizeof(double), cudaMemcpyDeviceToHost)); 

    //cudaMemcpy(calc, dev_calc, (n+1) * sizeof(double), cudaMemcpyDeviceToHost); 

    for (int i = 1; i <= n; i++) { 
     cout << "\nCALC cout = " << cuCreal(calc[i]) << ", " << cuCimag(calc[i]); 
    } 

    return 0; 
} 

Окончательный ответ неверен, и я также определил несколько других мест, где я не получаю ожидаемых значений.

1) Я ожидал сложный двойной массив (0, 5i) для всех элементов «j» после следующей строки кода. Тем не менее, я получаю все 0s. Почему это?

j[i] = make_cuDoubleComplex(0, 5); 

2) Почему я не могу напечатать свой массив с помощью cout? Строка кода, показанная ниже, дает следующую ошибку: ни один оператор «< <» не соответствует этим операндам. Как я могу исправить это, не используя printf?

cout << "\nJ = " << j[i]; 

3) Функция ГПУ 'FUNC', который должен выдавать массив (0, 10i) в качестве окончательного ответа дает случайные значения, такие как эти:

CALC = -1.#QNAN0 
CALC = -1.#QNAN0 
CALC = -9255963134931783100000000...000.. etc 
CALC = -9255963134931783100000000...000.. etc 

4) для данного фактическое исследование, комплексный массив «j» будет задан в формате сложного (двойного), а не cuDoubleComplex. Могу ли я выполнять аналогичные операции с массивом «j» сложных (двойных) с помощью функции «func»? Если нет, каковы мои варианты?

Я думаю, что я хорошо объяснил себя, но не стесняйтесь задавать любые последующие вопросы. Новое на C++, а также CUDA так неплохо: D

Answer 1

При написании кода CUDA, особенно, когда вы учитесь или с трудом (что-то не работает так, как вы ожидаете), вы должны всегда делать cuda error checking на всех CUDA Вызовы API и вызовы ядра.

Я не думаю, что на вашем коде есть функциональные ошибки CUDA (хорошая работа!), Но стоит отметить.

Большинство ваших вопросов связано с тем, что вы не печатаете тип cuDoubleComplex должным образом. В ваших заявлениях printf указывается параметр формата float (%f), но вы не передаете значение float (вы передаете значение cuDoubleComplex). Это не сработает, и printf будет вести себя странно, когда вы это сделаете, без указания каких-либо ошибок.

Вместо этого, попробовать что-то вроде этого:

printf("\nJ = %f, %f", cuCreal(j[i]), cuCimag(j[i])); 

Этих функции (cuCreal и cuCimag) возвращают действительные и мнимые части cuComplex чисел, и они возвращают их в качестве соответствующего типа, float или double, и в в этом случае неявный листинг от double до float подходит для того, что вы делаете, и может обрабатываться printf (хотя это не очень хорошая практика программирования - вместо этого используйте правильный спецификатор формата printf для значения double).

Если вы сделаете это изменение для обоих своих заявлений printf, я думаю, вы получите ожидаемые результаты - по крайней мере, я сделал, когда я запустил ваш код. Если вы по-прежнему получаете мусор, то ваш CUDA-графический процессор может работать некорректно, и вот где эта проверка ошибок CUDA, о которой я упоминал, поможет вам понять, в чем проблема.

Что касается ваших вопросов относительно cout, ответ примерно эквивалентен моему объяснению того, что происходит с printf. cout не понимает тип cuDoubleComplex и поэтому выдает ошибку. Если вы хотите исправить это, не используя printf, конвертируйте cuDoubleComplex в отдельные реальные и мнимые части, представленные float или double, используя функции преобразования, указанные в приведенном выше описании printf.

Что касается вашего последнего вопроса, то не должно быть сложно преобразовать ваши данные complex в cuDoubleComplex. Напишите функцию преобразования, чтобы сделать это на основе утилит, которые у вас есть в cuComplex.h. Есть обратные пути вокруг этого, но они не являются хорошей практикой программирования.

EDIT: В ответ на последующие вопросы было еще две ошибки в опубликованном коде.

1. dev_j и dev_calc имеют тип cuDoubleComplex но вы делали cudaMalloc и cudaMemcpy от этих величин, как если бы они были размером double. В следующем коде я изменил эти sizeof(double) записей на sizeof(cuDoubleComplex).
2. Ваша индексация вообще была немного странной для C и C++. Обычно индексы начинаются с нуля. У вас была проблема с индексацией, когда последний элемент не был правильно рассчитан. Я изменил всю индексацию на нулевую.

Вот модификация кода, который работает для меня:

//#include <complex> // not necessary for this code 
#include <iostream> 
#include <cmath> 
//#include "cuda.h" // not necessary when compiling with nvcc 
#include "math.h" 
#include "cuComplex.h" 

#define n 5 

using namespace std; 

#define gpuErrchk(ans) { gpuAssert((ans), __FILE__, __LINE__); } 
inline void gpuAssert(cudaError_t code, 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); 
    } 
} 

__global__ void func(double *s, cuDoubleComplex *j, cuDoubleComplex *calc) { 

    int tid = blockIdx.x; 

    calc[tid] = cuCmul(j[tid], make_cuDoubleComplex(*s, 0)); 

} 

int main(void) { 


    cuDoubleComplex calc[n+1], *dev_j, *dev_calc; 
    double *dev_s, s[n] = { 2.0, 2.0, 2.0, 2.0, 2.0 }; 
    //complex<double> j[n+1] 
    cuDoubleComplex j[n]; 

    for (int i = 0; i < n; i++) { 
     j[i] = make_cuDoubleComplex(0, 5); 
     cout << "\nJ cout = " << cuCreal(j[i]) << ", " << cuCimag(j[i]); 
    } 

    // allocate the memory on the GPU 
    cudaMalloc((void**)&dev_s, (n) * sizeof(double)); 
    cudaMalloc((void**)&dev_j, (n) * sizeof(cuDoubleComplex)); 
    cudaMalloc((void**)&dev_calc, (n) * sizeof(cuDoubleComplex)); 

    cudaMemcpy(dev_s, s, (n) * sizeof(double), cudaMemcpyHostToDevice); 
    cudaMemcpy(dev_j, j, (n) * sizeof(cuDoubleComplex), cudaMemcpyHostToDevice); 

    func<<<n,1>>>(dev_s, dev_j, dev_calc); 
    //kernel<<<1,1>>>(a_d); 
    gpuErrchk(cudaPeekAtLastError()); 
    gpuErrchk(cudaMemcpy(calc, dev_calc, (n) * sizeof(cuDoubleComplex), cudaMemcpyDeviceToHost)); 

    //cudaMemcpy(calc, dev_calc, (n+1) * sizeof(double), cudaMemcpyDeviceToHost); 

    for (int i = 0; i < n; i++) { 
     cout << "\nCALC cout = " << cuCreal(calc[i]) << ", " << cuCimag(calc[i]); 
    } 

    return 0; 
}