Массив комплексных чисел в PyopenCL

Question

Я работаю над новой проблемой с PyopenCl, в которой мне приходится иметь дело со сложными числами. Точнее, было бы очень удобно использовать массив 2D numpy со сложными числами внутри. Что-то вроде: np_array [np_array [C_number, C_number, ..], np_array [C_number, C_number, ..], ...]

Тогда для результатов я должен был бы простой 1D Numpy массив комплексных чисел ,

Я также заметил, что pyopencl видит комплексное число numpy как float2, для которого я использую float16 для массива данных, так как у меня есть около 8 чисел, с которыми нужно иметь дело.

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

Код, который я использую, следующий.

import pyopencl as cl 
import numpy as np 

ctx = cl.create_some_context() 
queue = cl.CommandQueue(ctx) 
MF = cl.mem_flags 

M = 3 

zero = np.complex64(0.0) 

X1_h = np.array([1 + 1j*2, 2 + 1j*3, 3 + 1j*4]).astype(np.complex64) 
X2_h = np.array([1 + 1j*2, 2 + 1j*3, 3 + 1j*4]).astype(np.complex64) 
X3_h = np.array([1 + 1j*2, 2 + 1j*3, 3 + 1j*4]).astype(np.complex64) 
Y1_h = np.array([4 + 1j*5, 5 + 1j*6, 6 + 1j*7]).astype(np.complex64) 
Y2_h = np.array([4 + 1j*5, 5 + 1j*6, 6 + 1j*7]).astype(np.complex64) 
Y3_h = np.array([4 + 1j*5, 5 + 1j*6, 6 + 1j*7]).astype(np.complex64) 
aux_h = np.complex64(1 + 1j*1) 
RES_h = np.empty_like(X1_h) 

dados_h = [] 
for i in range(3): 
    dados_h.append(np.array([X1_h[i], X2_h[i], X3_h[i], Y1_h[i], Y2_h[i], Y3_h[i]]).astype(np.complex64)) 
dados_h = np.array(dados_h).astype(np.complex64) 

print dados_h 

aux_d = cl.Buffer(ctx, MF.READ_WRITE | MF.COPY_HOST_PTR, hostbuf=aux_h) 
dados_d = cl.Buffer(ctx, MF.READ_WRITE | MF.COPY_HOST_PTR, hostbuf=dados_h) 
RES_d = cl.Buffer(ctx, MF.READ_WRITE | MF.COPY_HOST_PTR, hostbuf = RES_h) 

Source = """ 
__kernel void soma(__global float2 *aux, __global float16 *dados, __global float2 *res){ 
const int gid_x = get_global_id(0); 
const int gid_y = get_global_id(1); 
res[gid_x].x = dados[gid_y].s0; 
res[gid_x].y = dados[gid_y].s1; 
} 
""" 
prg = cl.Program(ctx, Source).build() 

completeEvent = prg.soma(queue, (M,), None, aux_d, dados_d, RES_d) 
completeEvent.wait() 

cl.enqueue_copy(queue, RES_h, RES_d) 
print "GPU" 
print RES_h 

Выход я получаю следующее:

[[ 1.+2.j 1.+2.j 1.+2.j 4.+5.j 4.+5.j 4.+5.j] 
[ 2.+3.j 2.+3.j 2.+3.j 5.+6.j 5.+6.j 5.+6.j] 
[ 3.+4.j 3.+4.j 3.+4.j 6.+7.j 6.+7.j 6.+7.j]] 
GPU 
[ 1.+2.j 1.+2.j 1.+2.j] 

Мой ожидается выход:

[ 1.+2.j 2.+3.j 3.+4.j] 

Я не могу понять, как я получаю этот результат. Как я уже сказал, я считаю, что это связано с идентификаторами потоков, но я не могу понять это. Если я использую gid_x вместо gid_y на правой части красного [gid_x] = ... Я получаю следующее

[ 1.+2.j 2.+3.j 6.+7.j] 

Может кто-нибудь дать мне некоторое представление в то, что я делаю неправильно, пожалуйста?

Answer 1

Вы запускаете одномерное ядро, поэтому get_global_id(1) всегда будет возвращать 0. Это объясняет, почему ваше ядро ​​просто копирует первый элемент массива dados в каждый элемент вывода.

Использование float16 для представления одной «строки» вашего ввода работает только в том случае, если у вас на самом деле имеется 8 комплексных чисел в строке. В вашем примере у вас всего 6, поэтому вы не получаете правильных результатов при копировании с dados[gid_x].

Чтобы ваш код, чтобы иметь дело с произвольным размером строки, просто передать ширину в качестве параметра, а затем вычислить линейный индекс вручную:

__kernel void soma(__global float2 *aux, 
        __global float2 *dados, 
        __global float2 *res, 
        int rowWidth){ 
    const int gid_x = get_global_id(0); 
    res[gid_x] = dados[gid_x*rowWidth]; 
} 

, а затем передать строку ширины, как дополнительный аргумент при запуске ядра:

# Pass your actual row-width instead of 6! 
completeEvent = prg.soma(queue, (M,), None, aux_d, dados_d, RES_d, np.int32(6))