OpenCL 嵌套循环产生意想不到的结果

Question

我对 OpenCL/C99 有点陌生，但不明白为什么下面的两个内核会给出不同的结果。 X 已初始化为零，但需要在每个外部循环步骤“重新归零”，否则会获得不正确的结果（参见图）。请注意，我在这里没有调用任何并行性；据我所知，这完全是连续的：

kernels.cl的内容：

#pragma OPENCL EXTENSION cl_khr_fp64: enable
#define __CL_ENABLE_EXCEPTIONS

__kernel void nested_sum_succeeds(
    int L,
    __global __read_only double* X,
    __global double* Y
)
{
    for (int k=0; k<=L; k++) {
        Y[k] = 0;
        for (int j=0; j<=k; j++) {
            Y[k] += X[j];
        }
    }
}

__kernel void nested_sum_fails(
    int L,
    __global __read_only double* X,
    __global double* Y
)
{
    for (int k=0; k<=L; k++) {
        // Y[k] = 0;
        for (int j=0; j<=k; j++) {
            Y[k] += X[j];
        }
    }
}

script.py的内容：

import numpy as np
import pyopencl as cl
import pyopencl.array as cl_array
import matplotlib.pyplot as plt

with open("./kernels.cl") as fp:
    prog_str = fp.read()
ctx = cl.create_some_context()
queue = cl.CommandQueue(ctx)
prog = cl.Program(ctx, prog_str).build()

L = 1000
X = np.linspace(0, 10, L)
X_dev = cl_array.to_device(queue, X)
Y_succeeds_dev = cl_array.to_device(queue, np.zeros(shape=X.shape, dtype=np.float64))
Y_fails_dev = cl_array.to_device(queue, np.zeros(shape=X.shape, dtype=np.float64))


nested_sum_succeeds = prog.nested_sum_succeeds
nested_sum_succeeds.set_scalar_arg_dtypes([
    np.int64,
    None,
    None,
])

nested_sum_succeeds(
    queue,
    (len(X),),
    None,
    L,
    X_dev.data,
    Y_succeeds_dev.data,
)


nested_sum_fails = prog.nested_sum_fails
nested_sum_fails.set_scalar_arg_dtypes([
    np.int64,
    None,
    None,
])

nested_sum_fails(
    queue,
    (len(X),),
    None,
    L,
    X_dev.data,
    Y_fails_dev.data,
)

np.allclose(Y_succeeds_dev.get(), Y_fails_dev.get()) #False

plt.ion()
plt.plot(Y_succeeds_dev.get())
plt.plot(Y_fails_dev.get())

结果：

Answer 1

请注意，我在这里没有调用任何并行性；

是的，内核确实是并行运行的——它计划在第一维运行len(x) 工作项。将其更改为使用 1 个工作项进行处理后一切正常：

nested_sum_succeeds(
    queue,
    (1,),
    None,
    L,
    X_dev.data,
    Y_succeeds_dev.data,
)

nested_sum_fails(
    queue,
    (1,),
    None,
    L,
    X_dev.data,
    Y_fails_dev.data,
)

然后np.allclose(Y_succeeds_dev.get(), Y_fails_dev.get()) 返回 True。 您还可以从 nested_sum_succeeds 内核中删除该归零 Y[k] = 0;，因为它不需要。

此外，如果您想在其他设备上运行此内核，则需要进行一些小修复，因为并非所有编译器都会接受第一个内核参数的类型在内核 int 中，而在 python 中指定为 @ 987654327@，必须匹配内核中的内容，所以：

nested_sum_succeeds.set_scalar_arg_dtypes([
    np.int32,
    None,
    None,
])

nested_sum_fails.set_scalar_arg_dtypes([
    np.int32,
    None,
    None,
])

还有一件事适用于在其他设备上使用，我将删除 __read_only 访问限定符，它也不会在所有设备上编译。