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 个数字要处理。

为了完成基本操作，我构建了一个简单的程序。 我已经制定了构建首字母数组并将它们发送到内核的方法，但结果与我的预期不同。我想这与线程的 ID 有关，但我似乎无法弄清楚。

我使用的代码如下。

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]

我不明白我是如何得到这个结果的。如前所述，我认为这与线程 ID 有关，但我无法弄清楚。如果我在 red[gid_x] = ... 的右侧使用 gid_x 而不是 gid_y，我会得到以下结果

[ 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))