PyopenCL 中的复数数组
Array of complex numbers in PyopenCL
我一直在处理 PyopenCl 的一个新问题,其中我必须处理复数。更准确地说,使用内部包含复数的 2D numpy 数组会非常方便。
就像是:
np_array[np_array[C_number, C_number, ..], np_array[C_number, C_number, ..], ...]
然后为了得到结果,我需要一个简单的一维 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 相关,但我无法弄清楚。如果我在红色 [gid_x] 的右侧使用 gid_x 而不是 gid_y = ... 我会得到以下结果
[ 1.+2.j 2.+3.j 6.+7.j]
任何人都可以告诉我我做错了什么吗?
您正在启动一维内核,因此 get_global_id(1)
将始终 return 0
。这解释了为什么您的内核只是将 dados
数组的第一个元素复制到输出的每个元素中。
使用 float16
表示输入的一个 'row' 只有在每行实际上有 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))
我一直在处理 PyopenCl 的一个新问题,其中我必须处理复数。更准确地说,使用内部包含复数的 2D numpy 数组会非常方便。 就像是: np_array[np_array[C_number, C_number, ..], np_array[C_number, C_number, ..], ...]
然后为了得到结果,我需要一个简单的一维 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 相关,但我无法弄清楚。如果我在红色 [gid_x] 的右侧使用 gid_x 而不是 gid_y = ... 我会得到以下结果
[ 1.+2.j 2.+3.j 6.+7.j]
任何人都可以告诉我我做错了什么吗?
您正在启动一维内核,因此 get_global_id(1)
将始终 return 0
。这解释了为什么您的内核只是将 dados
数组的第一个元素复制到输出的每个元素中。
使用 float16
表示输入的一个 'row' 只有在每行实际上有 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))