在 pyopencl 中添加复杂的 3D 数组
Posted
技术标签:
【中文标题】在 pyopencl 中添加复杂的 3D 数组【英文标题】:Addition of complex 3D arrays in pyopencl 【发布时间】:2020-07-03 14:24:18 【问题描述】:我正在尝试在 python 库 pyopencl 的帮助下添加 2 个三维复杂数组。在 gpu 上获得的结果与使用 cpu 获得的参考结果不同。为什么 GPU 代码不提供复数的正确加法?我希望相应的变量 res 和 ref 彼此相等。我使用的代码:
import pyopencl as cl
import numpy as np
from numpy.fft import fftn, ifftn
from numpy.random import random, normal
import os
os.environ['PYOPENCL_COMPILER_OUTPUT'] = '1'
ctx = cl.create_some_context()
queue = cl.CommandQueue(ctx)
Lx = 100
Ly = 100
Lz = 1
L = Lx * Ly * Lz
const = np.zeros(4).astype(np.float32)
const[0] = Lx
const[1] = Ly
const[2] = Lz
const[3] = L
const_buf = cl.Buffer(ctx, cl.mem_flags.READ_ONLY | cl.mem_flags.COPY_HOST_PTR, hostbuf=const)
arr1 = np.random.rand(Lz, Ly, Lx).astype(np.float32)
arr2 = np.random.rand(Lz, Ly, Lx).astype(np.float32)
F1 = fftn(arr1)
F2 = fftn(arr2)
out = np.zeros((Lz,Ly,Lx)).astype(np.complex64)
out_buf = cl.Buffer(ctx, cl.mem_flags.WRITE_ONLY, out.nbytes)
do_it_code = """
#include <pyopencl-complex.h>
__kernel void doit(__global const float *constants,
__global const cfloat_t *IN1,__global const cfloat_t *IN2,
__global cfloat_t *out)
int z = get_global_id(0);
int y = get_global_id(1);
int x = get_global_id(2);
const int len = constants[3];
const int lx = constants[0];
const int ly = constants[1];
const int lz = constants[2];
const int pl = lx * ly;
int i = x + y * lx + z * pl;
if (x <= lx-1 && y <= ly-1 && z <= lz-1)
out[i] = cfloat_add(IN1[i],IN2[i]);
;
;
"""
bld_do_it = cl.Program(ctx, do_it_code).build()
def do_it(a,b):
a_buf = cl.Buffer(ctx, cl.mem_flags.READ_ONLY | cl.mem_flags.COPY_HOST_PTR, hostbuf=a.astype(np.complex64))
b_buf = cl.Buffer(ctx, cl.mem_flags.READ_ONLY | cl.mem_flags.COPY_HOST_PTR, hostbuf=b.astype(np.complex64))
launch = bld_do_it.doit(queue, a.shape, None, const_buf, a_buf, b_buf, out_buf)
launch.wait()
cl.enqueue_copy(queue, out, out_buf)
return out
ref=F1+F2
print(ref)
print("_________")
res=do_it(F1,F2)
print(res)
【问题讨论】:
结果的差异似乎是逗号后面的几位数字,这表明这可能是由于 CPU 和加速器上的不同 FPU 单元造成的。查看更多here 不,区别不仅仅是逗号后面的几个数字。 是的,你是对的,还有更多——看我的回答。 【参考方案1】:找出问题所在通常最好处理更小的数据集。让我们设置Lx=Ly=2
并检查我们得到了什么:
[[[ 3.40425836+0.j -2.09837677+0.j]
[-0.27708016+0.j 0.60454108+0.j]]]
_________
[[[ 3.4042583 +0.j -0.27708015+0.j]
[-2.0983768 +0.j 0.60454106+0.j]]]
数字几乎相同,只是并非都在预期的位置。简单的转置可以解决这个问题:
res=np.transpose(res, (0,2,1))
然后结果与参考结果匹配:
[[[ 3.4042583 +0.j -2.0983768 +0.j]
[-0.27708015+0.j 0.60454106+0.j]]]
结果之间的差异仍然很小 - 在逗号后的几位数字之后 - 但这可以通过 CPU 和 GPU 之间使用的不同 FPU 单元来解释(在我的例子中)。可以阅读更多相关信息,例如 here。
这不应该需要转置 - Lx,Ly,Lz 的放置一定是错误的 - 我会让 OP 弄清楚这一点。
【讨论】:
谢谢你,这很有道理!【参考方案2】:我发现对 OP 的 GPU 代码的此更新解决了讨论线程中提到的转置错误 - 并在一个回复中作为“学生练习”留下。事实证明,至少对我来说,弄清楚这件事并不简单,所以我认为这可能有助于分享:
do_it_code = """
#include <pyopencl-complex.h>
__kernel void doit(__global const float *constants,
__global const cfloat_t *IN1,__global const cfloat_t *IN2,
__global cfloat_t *out)
int z = get_global_id(0);
int y = get_global_id(1);
int x = get_global_id(2);
const int len = constants[3];
const int lx = constants[0];
const int ly = constants[1];
const int lz = constants[2];
const int pl = lx * ly;
const int pl2 = lz * ly;
int i = x + y * lx + z * pl;
int i2 = x*pl2 + y*lz + z;
if (x <= lx-1 && y <= ly-1 && z <= lz-1)
out[i] = cfloat_add(IN1[i2],IN2[i2]);
;
;
"""
【讨论】:
正如目前所写,您的答案尚不清楚。请edit 添加其他详细信息,以帮助其他人了解这如何解决所提出的问题。你可以找到更多关于如何写好答案的信息in the help center。以上是关于在 pyopencl 中添加复杂的 3D 数组的主要内容,如果未能解决你的问题,请参考以下文章