在 pyopencl 中添加复杂的 3D 数组

Question

我正在尝试在 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)

Answer 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 弄清楚这一点。

Answer 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]);

};
};
"""