OpenCL：2-10x 运行 时间求和列而不是方形阵列的行？
OpenCL: 2-10x run time summing columns rather that rows of square array?
我刚刚开始使用 OpenCL，所以我确信我可以做很多事情来改进这段代码，但有一件事对我来说特别突出：如果我对列而不是行求和 (基本上是连续的与跨步的，因为所有缓冲区都是线性的）在二维数据数组中，我得到不同的 运行 次，乘以 2 到 10 倍的任何地方。奇怪的是，连续访问似乎变慢了。
我正在使用 PyOpenCL 进行测试。
这是两个感兴趣的内核（reduce 和 reduce2），另一个正在生成为它们提供数据的内核（forcesCL）：
kernel void forcesCL(global float4 *chrgs, global float4 *chrgs2,float k, global float4 *frcs)
{
int i=get_global_id(0);
int j=get_global_id(1);
int idx=i+get_global_size(0)*j;
float3 c=chrgs[i].xyz-chrgs2[j].xyz;
float l=length(c);
frcs[idx].xyz= (l==0 ? 0
: ((chrgs[i].w*chrgs2[j].w)/(k*pown(l,2)))*normalize(c));
frcs[idx].w=0;
}
kernel void reduce(global float4 *frcs,ulong k,global float4 *result)
{
ulong gi=get_global_id(0);
ulong gs=get_global_size(0);
float3 tmp=0;
for(ulong i=0;i<k;i++)
tmp+=frcs[gi+i*gs].xyz;
result[gi].xyz=tmp;
}
kernel void reduce2(global float4 *frcs,ulong k,global float4 *result)
{
ulong gi=get_global_id(0);
ulong gs=get_global_size(0);
float3 tmp=0;
for(ulong i=0;i<k;i++)
tmp+=frcs[gi*gs+i].xyz;
result[gi].xyz=tmp;
}
这里感兴趣的是 reduce 内核。 forcesCL 内核仅估计两个电荷之间的洛伦兹力，其中每个电荷的位置在 float4 的 xyz 分量中编码，而电荷在 w 分量中。物理学不重要，OpenCL 只是一个玩具
除了展示构建步骤外，除非有人问我，否则我不会完成 PyOpenCL 设置：
program=cl.Program(context,'\n'.join((src_forcesCL,src_reduce,src_reduce2))).build()
然后我设置具有随机位置和基本电荷的阵列：
a=np.random.rand(10000,4).astype(np.float32)
a[:,3]=np.float32(q)
b=np.random.rand(10000,4).astype(np.float32)
b[:,3]=np.float32(q)
设置从头 space 和 return 值的分配：
c=np.empty((10000,10000,4),dtype=np.float32)
cc=cl.Buffer(context,cl.mem_flags.READ_WRITE,c.nbytes)
r=np.empty((10000,4),dtype=np.float32)
rr=cl.Buffer(context,cl.mem_flags.WRITE_ONLY,r.nbytes)
s=np.empty((10000,4),dtype=np.float32)
ss=cl.Buffer(context,cl.mem_flags.WRITE_ONLY,s.nbytes)
然后我尝试 运行 这两种方式中的每一种——一次使用 reduce()，一次使用 reduce2()。唯一的区别应该是我在哪个轴上求和：
%%timeit
aa=cl.Buffer(context,cl.mem_flags.READ_ONLY|cl.mem_flags.COPY_HOST_PTR,hostbuf=a)
bb=cl.Buffer(context,cl.mem_flags.READ_ONLY|cl.mem_flags.COPY_HOST_PTR,hostbuf=b)
evt1=program.forcesCL(queue,c.shape[0:2],None,aa,bb,k,cc)
evt2=program.reduce(queue,r.shape[0:1],None,cc,np.uint32(b.shape[0:1]),rr,wait_for=[evt1])
evt4=cl.enqueue_copy(queue,r,rr,wait_for=[evt2],is_blocking=True)
请注意，我已将参数交换为 forcesCL，因此我可以根据第一种方法检查结果：
%%timeit
aa=cl.Buffer(context,cl.mem_flags.READ_ONLY|cl.mem_flags.COPY_HOST_PTR,hostbuf=a)
bb=cl.Buffer(context,cl.mem_flags.READ_ONLY|cl.mem_flags.COPY_HOST_PTR,hostbuf=b)
evt1=program.forcesCL(queue,c.shape[0:2],None,bb,aa,k,cc)
evt2=program.reduce2(queue,s.shape[0:1],None,cc,np.uint32(a.shape[0:1]),ss,wait_for=[evt1])
evt4=cl.enqueue_copy(queue,s,ss,wait_for=[evt2],is_blocking=True)
使用 reduce() 内核的版本给我大约 140 毫秒的时间，使用 reduce2() 内核的版本给我大约 360 毫秒的时间。 returned 的值是相同的，只是符号更改，因为它们是以相反的顺序被减去的。
如果我执行一次 forcesCL 步骤，然后 运行 执行两个 reduce 内核，差异会更加明显——大约 30 毫秒与 250 毫秒。
我没想到会有什么不同，但如果我是的话，我会期望连续访问的性能更好，而不是更差。
谁能告诉我这里发生了什么？
谢谢！
这是合并的一个明显例子。这与索引的使用方式无关（在行或列中），而是在 HW 中如何访问内存。只需逐步了解实际访问的执行方式和顺序即可。
让我们好好分析一下：
假设工作项被划分为大小为 N.
的本地块
第一种情况：
WI_0 will read 0, Gs, 2Gs, 3Gs, .... (k-1)Gs
WI_1 will read 1, Gs+1, 2Gs+1, 3Gs+1, .... (k-1)Gs+1
...
由于每个WI都是运行并行的，所以它们的内存访问是同时发生的。因此，请求内存控制器：
First iteration: 0, 1, 2, 3 ... N-1 -> Groups into few memory access
Second iteration: Gs, Gs+1, Gs+2, ... Gs+N-1 -> Groups into few memory access
...
在那种情况下，在每次迭代中，内存控制器将所有 N WI 请求分组到一个大的 256 位 reads/writes 到全局。不需要缓存，因为处理完数据后可以丢弃。
第二种情况：
WI_0 will read 0, 1, 2, 3, .... (k-1)
WI_1 will read Gs, Gs+1, Gs+2, Gs+3, .... Gs+(k-1)
...
请求内存控制器：
First iteration: 0, Gs, 2Gs, 3Gs -> Scattered read, no grouping
Second iteration: 1, Gs+1, 2Gs+1, 3Gs+1 ->Scattered read, no grouping
...
在这种情况下，内存控制器未在正确的模式下工作。如果缓存内存是无限的，它会起作用，但事实并非如此。由于有时请求的工作项间内存相同（由于循环的 k 大小），它可以缓存一些读取，但不是全部。
当您减小 k 的值时，您就减少了可能的缓存重用量 I。这导致列和行访问模式之间的差异更大。
我刚刚开始使用 OpenCL，所以我确信我可以做很多事情来改进这段代码，但有一件事对我来说特别突出：如果我对列而不是行求和 (基本上是连续的与跨步的，因为所有缓冲区都是线性的）在二维数据数组中，我得到不同的 运行 次，乘以 2 到 10 倍的任何地方。奇怪的是，连续访问似乎变慢了。
我正在使用 PyOpenCL 进行测试。
这是两个感兴趣的内核（reduce 和 reduce2），另一个正在生成为它们提供数据的内核（forcesCL）：
kernel void forcesCL(global float4 *chrgs, global float4 *chrgs2,float k, global float4 *frcs)
{
int i=get_global_id(0);
int j=get_global_id(1);
int idx=i+get_global_size(0)*j;
float3 c=chrgs[i].xyz-chrgs2[j].xyz;
float l=length(c);
frcs[idx].xyz= (l==0 ? 0
: ((chrgs[i].w*chrgs2[j].w)/(k*pown(l,2)))*normalize(c));
frcs[idx].w=0;
}
kernel void reduce(global float4 *frcs,ulong k,global float4 *result)
{
ulong gi=get_global_id(0);
ulong gs=get_global_size(0);
float3 tmp=0;
for(ulong i=0;i<k;i++)
tmp+=frcs[gi+i*gs].xyz;
result[gi].xyz=tmp;
}
kernel void reduce2(global float4 *frcs,ulong k,global float4 *result)
{
ulong gi=get_global_id(0);
ulong gs=get_global_size(0);
float3 tmp=0;
for(ulong i=0;i<k;i++)
tmp+=frcs[gi*gs+i].xyz;
result[gi].xyz=tmp;
}
这里感兴趣的是 reduce 内核。 forcesCL 内核仅估计两个电荷之间的洛伦兹力，其中每个电荷的位置在 float4 的 xyz 分量中编码，而电荷在 w 分量中。物理学不重要，OpenCL 只是一个玩具
除了展示构建步骤外，除非有人问我，否则我不会完成 PyOpenCL 设置：
program=cl.Program(context,'\n'.join((src_forcesCL,src_reduce,src_reduce2))).build()
然后我设置具有随机位置和基本电荷的阵列：
a=np.random.rand(10000,4).astype(np.float32)
a[:,3]=np.float32(q)
b=np.random.rand(10000,4).astype(np.float32)
b[:,3]=np.float32(q)
设置从头 space 和 return 值的分配：
c=np.empty((10000,10000,4),dtype=np.float32)
cc=cl.Buffer(context,cl.mem_flags.READ_WRITE,c.nbytes)
r=np.empty((10000,4),dtype=np.float32)
rr=cl.Buffer(context,cl.mem_flags.WRITE_ONLY,r.nbytes)
s=np.empty((10000,4),dtype=np.float32)
ss=cl.Buffer(context,cl.mem_flags.WRITE_ONLY,s.nbytes)
然后我尝试 运行 这两种方式中的每一种——一次使用 reduce()，一次使用 reduce2()。唯一的区别应该是我在哪个轴上求和：
%%timeit
aa=cl.Buffer(context,cl.mem_flags.READ_ONLY|cl.mem_flags.COPY_HOST_PTR,hostbuf=a)
bb=cl.Buffer(context,cl.mem_flags.READ_ONLY|cl.mem_flags.COPY_HOST_PTR,hostbuf=b)
evt1=program.forcesCL(queue,c.shape[0:2],None,aa,bb,k,cc)
evt2=program.reduce(queue,r.shape[0:1],None,cc,np.uint32(b.shape[0:1]),rr,wait_for=[evt1])
evt4=cl.enqueue_copy(queue,r,rr,wait_for=[evt2],is_blocking=True)
请注意，我已将参数交换为 forcesCL，因此我可以根据第一种方法检查结果：
%%timeit
aa=cl.Buffer(context,cl.mem_flags.READ_ONLY|cl.mem_flags.COPY_HOST_PTR,hostbuf=a)
bb=cl.Buffer(context,cl.mem_flags.READ_ONLY|cl.mem_flags.COPY_HOST_PTR,hostbuf=b)
evt1=program.forcesCL(queue,c.shape[0:2],None,bb,aa,k,cc)
evt2=program.reduce2(queue,s.shape[0:1],None,cc,np.uint32(a.shape[0:1]),ss,wait_for=[evt1])
evt4=cl.enqueue_copy(queue,s,ss,wait_for=[evt2],is_blocking=True)
使用 reduce() 内核的版本给我大约 140 毫秒的时间，使用 reduce2() 内核的版本给我大约 360 毫秒的时间。 returned 的值是相同的，只是符号更改，因为它们是以相反的顺序被减去的。
如果我执行一次 forcesCL 步骤，然后 运行 执行两个 reduce 内核，差异会更加明显——大约 30 毫秒与 250 毫秒。
我没想到会有什么不同，但如果我是的话，我会期望连续访问的性能更好，而不是更差。
谁能告诉我这里发生了什么？
谢谢！
这是合并的一个明显例子。这与索引的使用方式无关（在行或列中），而是在 HW 中如何访问内存。只需逐步了解实际访问的执行方式和顺序即可。
让我们好好分析一下：
假设工作项被划分为大小为 N.
的本地块第一种情况：
WI_0 will read 0, Gs, 2Gs, 3Gs, .... (k-1)Gs
WI_1 will read 1, Gs+1, 2Gs+1, 3Gs+1, .... (k-1)Gs+1
...
由于每个WI都是运行并行的，所以它们的内存访问是同时发生的。因此，请求内存控制器：
First iteration: 0, 1, 2, 3 ... N-1 -> Groups into few memory access
Second iteration: Gs, Gs+1, Gs+2, ... Gs+N-1 -> Groups into few memory access
...
在那种情况下，在每次迭代中，内存控制器将所有 N WI 请求分组到一个大的 256 位 reads/writes 到全局。不需要缓存，因为处理完数据后可以丢弃。
第二种情况：
WI_0 will read 0, 1, 2, 3, .... (k-1)
WI_1 will read Gs, Gs+1, Gs+2, Gs+3, .... Gs+(k-1)
...
请求内存控制器：
First iteration: 0, Gs, 2Gs, 3Gs -> Scattered read, no grouping
Second iteration: 1, Gs+1, 2Gs+1, 3Gs+1 ->Scattered read, no grouping
...
在这种情况下，内存控制器未在正确的模式下工作。如果缓存内存是无限的，它会起作用，但事实并非如此。由于有时请求的工作项间内存相同（由于循环的 k 大小），它可以缓存一些读取，但不是全部。
当您减小 k 的值时，您就减少了可能的缓存重用量 I。这导致列和行访问模式之间的差异更大。