如何解释这个 GPU 计算？
How to explain this GPU computation?
我正在使用 OpenCL 在我的 GPU 上实现二维图像处理内核。我从我的 GPU 得到了非常令人费解的结果。该代码使用 2X2 模板并计算模板中每个输入样本的平均值，并将计算出的平均值添加到位于模板内的输出图像中的每个样本。
这是 CPU 的代码：
for(int i2 = 1; i2 < n1; ++i2) {
for(int i1 = 1; i1 < n2; ++i1) {
r00 = h_r[i2 ][i1 ];
r01 = h_r[i2 ][i1-1];
r10 = h_r[i2-1][i1 ];
r11 = h_r[i2-1][i1-1];
rs = 0.25f*(r00+r01+r10+r11);
s[i2 ][i1 ] += rs;
s[i2 ][i1-1] += rs;
s[i2-1][i1 ] += rs;
s[i2-1][i1-1] += rs;
}
}
使用
2 2 2 2
2 2 2 2
2 2 2 2
2 2 2 2
作为输入图像，我在应用这个内核后得到以下输出图像：
2 4 4 2
4 8 8 4
4 8 8 4
2 4 4 2
在我的 OpenCL 实现中，我有以下内核：
_kernel void soSmoothingNew(__global const float* restrict d_r,
__global float* restrict d_s,
int n1,
int n2,)
{
int g1 = get_global_id(0);
int g0 = get_global_id(1);
int i1 = g1+1;
int i2 = g0+1;
if (i1 >= n2) return;
if (i2 >= n1) return;
float r00, r01, r10, r11, rs;
r01 = d_r[i2*n2+(i1-1)];
r10 = d_r[(i2-1)*n2 + i1];
r11 = d_r[(i2-1)*n2 + (i1-1)];
rs = 0.25f*(r00+r01+r10+r11);
d_s[i2*n2 + i1] += rs;
d_s[i2*n2 + (i1-1)] += rs;
d_s[(i2-1)*n2 + i1] += rs;
d_s[(i2-1)*n2 + (i1-1)] += rs;
}
结果输出为：
2 2 2 2
4 6 6 4
4 6 6 4
2 4 4 2
我正在使用以下代码从主机执行内核：
size_t local_group_size[2] = {4,4};
size_t global_group_size_block[2] = {ceil((n1/local_group_size[0]) + 1) * local_group_size[0],
ceil((n2/local_group_size[1]) + 1) * local_group_size[1]};
err = clEnqueueNDRangeKernel(queue, kernel1, 2, NULL, global_group_size_block, local_group_size, 0, NULL, NULL);
为了简洁起见，我省略了 clSetKernelArg、clCreateBuffer 和其他 OpenCL 调用。请注意，我还有另一个内核，它在执行该内核之前将 GPU 上的输出 d_s 数组清零。
我很难理解 GPU 上的线程是如何运行以达到这个结果的。对此有任何见解将不胜感激。
如void_ptr所述，您的问题绝对是写回输出时的竞争条件。解决这个问题的一个简单方法是让每个工作项负责完全计算其输出像素。你的算法也可以在这个过程中得到简化。
每个像素位于四个 2x2 像素框内，构成图像的 3x3 区域。
O O O
O X O X is the one we want to compute with a work item.
O O O
2x2 区域可以称为 As、Bs、Cs 和 Ds：
A A O O B B O O O O O O
A A O O B B C C O O D D
O O O O O O C C O O D D
您可能已经注意到，在使用您的原始算法时，X 像素是 4 个独立平均值的一部分。生成的像素恰好包含 1.0f * X。每个周围的像素都可以以类似的方式进行加权，并使用此掩码添加到最终输出值：
0.25 0.50 0.25
0.50 1.00 0.50
0.25 0.50 0.25
这些值有助于得出下面的内核。
_kernel void soSmoothing(__global const float* restrict d_r, __global float* restrict d_s, int n1, int n2)
{
//(idX,idY) represents the position to compute and write to d_s
int idX = get_global_id(0);
int idY = get_global_id(1);
//using zero-base indices
if(idX >= n1) return;
if(idY >= n2) return;
float outValue = 0.0f;
if(idX > 0){
outValue += 0.50f * d_r[idX-1 + idY*n1] + 0.25 * d_r[idX + idY*n1];
if (idY > 0){
outValue += 0.25f * d_r[idX-1 + (idY-1)*n1];
}
if (idY < (n2-1)){
outValue += 0.25f * d_r[idX-1 + (idY+1)*n1];
}
}
if(idX < (n1-1)){
outValue += 0.50f * d_r[idX+1 + idY*n1] + 0.25 * d_r[idX + idY*n1];
if (idY > 0){
outValue += 0.25f * d_r[idX+1 + (idY-1)*n1];
}
if (idY < (n2-1)){
outValue += 0.25f * d_r[idX+1 + (idY+1)*n1];
}
}
if (idY > 0){
outValue += 0.50f * d_r[idX + (idY-1)*n1] + 0.25 * d_r[idX + idY*n1];
}
if (idY < (n2-1)){
outValue += 0.50f * d_r[idX + (idY+1)*n1] + 0.25 * d_r[idX + idY*n1];
}
d_s[idX + idY*n1] = outValue;
}
这绝不是计算输出的最快方法，但它会为您提供正确的结果。通过仅向输出的每个元素写入一次来解决写入全局内存的竞争条件问题。
还有很大的改进空间
- 它不是最优的一个重要原因是每个像素平均将被读取 9 次——它周围的每个工作项一次，它自己的工作项一次。如果您遇到这个问题，我可以 post 更新来解决这个问题。
- 使运行更快的第二种方法是避免条件检查。为此，您可以在输入周围添加一个 1 像素的黑色边框，并在 x 和 y 轴上将坐标偏移调整为 +1。我在下面以'+n1+1'的形式进行了偏移。
可以避免检查，因为您知道坐标将落在输入图像的范围内。内核的其余部分保持不变。
float outValue = d_r[idX + idY*n1 + n1+1];
outValue += 0.50f * d_r[idX-1 + idY*n1 + n1+1];
outValue += 0.25f * d_r[idX-1 + (idY-1)*n1 + n1+1];
outValue += 0.25f * d_r[idX-1 + (idY+1)*n1 + n1+1];
outValue += 0.50f * d_r[idX+1 + idY*n1 + n1+1];
outValue += 0.25f * d_r[idX+1 + (idY-1)*n1 + n1+1];
outValue += 0.25f * d_r[idX+1 + (idY+1)*n1 + n1+1];
outValue += 0.50f * d_r[idX + (idY-1)*n1 + n1+1];
outValue += 0.50f * d_r[idX + (idY+1)*n1 + n1+1];
这是您询问的 "rows and columns" 方法。这个想法是同时处理输出图像的非重叠区域。
我使用的 16x16 像素示例可以扩展到任何图像尺寸。每个黑框代表一个像素。不同深浅的蓝色仅用于区分要计算的工作项区域。
如果将图像分成两个像素宽的列，则可以在单独的线程中处理每一列。这可以防止列之间的全局写入冲突。对奇数列重复 - 即 dx=1 以计算与列重叠的 2x2 分组，如第二个图像所示。这种方法只是防止了水平方向的写冲突，但是还是需要考虑到行数。
为了避免列内的写入冲突，您需要进一步将列分成行，如上图所示。让我们让一个工作项计算单个 2x2 输出以使内核尽可能简单。行也是分两个阶段计算的——dy=0 和 dy=1。总的来说，这是一个四步算法，其中每一步都可以称为 'embarrassingly parallel'。这些步骤可以通过它们相对于原点的位置偏移来引用：(dx,dy) = A(0,0)、B(1,0)、C(1,0) 和 D(1,1)。请注意，这些组不一定使用相同数量的工作项，具体取决于最终 row/column 是否足以进行计算。
显示内核调用期间两个示例工作项负责的插图。 "rows and columns" 基本上变成了“4 个棋盘”。下面是内核，它将计算 2x2 像素组的平均值并添加全局缓冲区。只要A、B、C、D组各自独立执行，就不需要检查全局写冲突。
_kernel void soSmoothing(__global const float* restrict inData, __global float* restrict outData, int width, int height, int dx, int dy)
{
//(idX,idY) represents the position to compute and write to outData
int idX = get_global_id(0);
int idY = get_global_id(1);
//(pixelX, pixelY) is the top-left corner of the square to compute
int pixelX = idX *2 +dx;
int pixelY = idY *2 +dy;
if(pixelX > (width-2)) return;
if(pixelY > (height-2)) return;
//do the math and add to the corresponding addresses in outData
int topLeftAddress = pixelX + pixelY * width;
float avg = 0.25f * (inData[topLeftAddress] + inData[topLeftAddress +1] + inData[topLeftAddress + width] + inData[topLeftAddress + width +1]);
outData[topLeftAddress] = outData[topLeftAddress] + avg;
outData[topLeftAddress +1] = outData[topLeftAddress +1] + avg;
outData[topLeftAddress +width] = outData[topLeftAddress +width] + avg;
outData[topLeftAddress +width +1] = outData[topLeftAddress +width +1] + avg;
}
主机代码步骤：
1) 将命令队列设置为阻塞
2) 设置缓冲区：1x 输入，1x 输出
3) 将输入缓冲区复制到设备，在设备上将输出初始化为 0
4a) 使用 dx=0, dy=0
使内核入队
4b) 使用 dx=1, dy=0
使内核入队
4c) 使用 dx=0, dy=1
使内核入队
4d) 使用 dx=1, dy=1
使内核入队
5) 将输出缓冲区读回主机
请注意，在所有内核执行完成之前，不需要将输出复制到主机。
优点：
- 所有 4 个内核的相同输入 运行 意味着您可以非常轻松地 运行 在最多 4 个独立的设备上执行此操作，并在主机上汇总它们的输出。也许尝试 运行 在 cpu 上安装 1 个内核，在 gpu 上安装 3 个内核。您还可以分解内核以处理图像的较小区域，并根据需要让尽可能多的 opencl 设备参与。传输的开销最终会成为主要瓶颈，因此这可能不值得付出努力。
- 内核简单，易于调试。
- 内核是 'chained'，因此它们会添加到之前的输出中。可以是 运行 任意顺序。
- 内核（+主机代码）可以轻松扩展以支持 3x3 或其他单元尺寸。
缺点：
- 输入图像被读取 4x
- 输出像素写入 4x
- 每个工作项一次只计算 4 个像素。计算更大的区域可以更快，但代码复杂度更高。
给你。请让我知道这对您来说效果如何，以及我是否需要进行任何更正。
我正在使用 OpenCL 在我的 GPU 上实现二维图像处理内核。我从我的 GPU 得到了非常令人费解的结果。该代码使用 2X2 模板并计算模板中每个输入样本的平均值，并将计算出的平均值添加到位于模板内的输出图像中的每个样本。
这是 CPU 的代码：
for(int i2 = 1; i2 < n1; ++i2) {
for(int i1 = 1; i1 < n2; ++i1) {
r00 = h_r[i2 ][i1 ];
r01 = h_r[i2 ][i1-1];
r10 = h_r[i2-1][i1 ];
r11 = h_r[i2-1][i1-1];
rs = 0.25f*(r00+r01+r10+r11);
s[i2 ][i1 ] += rs;
s[i2 ][i1-1] += rs;
s[i2-1][i1 ] += rs;
s[i2-1][i1-1] += rs;
}
}
使用
2 2 2 2
2 2 2 2
2 2 2 2
2 2 2 2
作为输入图像，我在应用这个内核后得到以下输出图像：
2 4 4 2
4 8 8 4
4 8 8 4
2 4 4 2
在我的 OpenCL 实现中，我有以下内核：
_kernel void soSmoothingNew(__global const float* restrict d_r,
__global float* restrict d_s,
int n1,
int n2,)
{
int g1 = get_global_id(0);
int g0 = get_global_id(1);
int i1 = g1+1;
int i2 = g0+1;
if (i1 >= n2) return;
if (i2 >= n1) return;
float r00, r01, r10, r11, rs;
r01 = d_r[i2*n2+(i1-1)];
r10 = d_r[(i2-1)*n2 + i1];
r11 = d_r[(i2-1)*n2 + (i1-1)];
rs = 0.25f*(r00+r01+r10+r11);
d_s[i2*n2 + i1] += rs;
d_s[i2*n2 + (i1-1)] += rs;
d_s[(i2-1)*n2 + i1] += rs;
d_s[(i2-1)*n2 + (i1-1)] += rs;
}
结果输出为：
2 2 2 2
4 6 6 4
4 6 6 4
2 4 4 2
我正在使用以下代码从主机执行内核：
size_t local_group_size[2] = {4,4};
size_t global_group_size_block[2] = {ceil((n1/local_group_size[0]) + 1) * local_group_size[0],
ceil((n2/local_group_size[1]) + 1) * local_group_size[1]};
err = clEnqueueNDRangeKernel(queue, kernel1, 2, NULL, global_group_size_block, local_group_size, 0, NULL, NULL);
为了简洁起见，我省略了 clSetKernelArg、clCreateBuffer 和其他 OpenCL 调用。请注意，我还有另一个内核，它在执行该内核之前将 GPU 上的输出 d_s 数组清零。
我很难理解 GPU 上的线程是如何运行以达到这个结果的。对此有任何见解将不胜感激。
如void_ptr所述，您的问题绝对是写回输出时的竞争条件。解决这个问题的一个简单方法是让每个工作项负责完全计算其输出像素。你的算法也可以在这个过程中得到简化。
每个像素位于四个 2x2 像素框内，构成图像的 3x3 区域。
O O O
O X O X is the one we want to compute with a work item.
O O O
2x2 区域可以称为 As、Bs、Cs 和 Ds：
A A O O B B O O O O O O
A A O O B B C C O O D D
O O O O O O C C O O D D
您可能已经注意到，在使用您的原始算法时，X 像素是 4 个独立平均值的一部分。生成的像素恰好包含 1.0f * X。每个周围的像素都可以以类似的方式进行加权，并使用此掩码添加到最终输出值：
0.25 0.50 0.25
0.50 1.00 0.50
0.25 0.50 0.25
这些值有助于得出下面的内核。
_kernel void soSmoothing(__global const float* restrict d_r, __global float* restrict d_s, int n1, int n2)
{
//(idX,idY) represents the position to compute and write to d_s
int idX = get_global_id(0);
int idY = get_global_id(1);
//using zero-base indices
if(idX >= n1) return;
if(idY >= n2) return;
float outValue = 0.0f;
if(idX > 0){
outValue += 0.50f * d_r[idX-1 + idY*n1] + 0.25 * d_r[idX + idY*n1];
if (idY > 0){
outValue += 0.25f * d_r[idX-1 + (idY-1)*n1];
}
if (idY < (n2-1)){
outValue += 0.25f * d_r[idX-1 + (idY+1)*n1];
}
}
if(idX < (n1-1)){
outValue += 0.50f * d_r[idX+1 + idY*n1] + 0.25 * d_r[idX + idY*n1];
if (idY > 0){
outValue += 0.25f * d_r[idX+1 + (idY-1)*n1];
}
if (idY < (n2-1)){
outValue += 0.25f * d_r[idX+1 + (idY+1)*n1];
}
}
if (idY > 0){
outValue += 0.50f * d_r[idX + (idY-1)*n1] + 0.25 * d_r[idX + idY*n1];
}
if (idY < (n2-1)){
outValue += 0.50f * d_r[idX + (idY+1)*n1] + 0.25 * d_r[idX + idY*n1];
}
d_s[idX + idY*n1] = outValue;
}
这绝不是计算输出的最快方法，但它会为您提供正确的结果。通过仅向输出的每个元素写入一次来解决写入全局内存的竞争条件问题。
还有很大的改进空间
- 它不是最优的一个重要原因是每个像素平均将被读取 9 次——它周围的每个工作项一次，它自己的工作项一次。如果您遇到这个问题，我可以 post 更新来解决这个问题。
- 使运行更快的第二种方法是避免条件检查。为此，您可以在输入周围添加一个 1 像素的黑色边框，并在 x 和 y 轴上将坐标偏移调整为 +1。我在下面以'+n1+1'的形式进行了偏移。
可以避免检查，因为您知道坐标将落在输入图像的范围内。内核的其余部分保持不变。
float outValue = d_r[idX + idY*n1 + n1+1];
outValue += 0.50f * d_r[idX-1 + idY*n1 + n1+1];
outValue += 0.25f * d_r[idX-1 + (idY-1)*n1 + n1+1];
outValue += 0.25f * d_r[idX-1 + (idY+1)*n1 + n1+1];
outValue += 0.50f * d_r[idX+1 + idY*n1 + n1+1];
outValue += 0.25f * d_r[idX+1 + (idY-1)*n1 + n1+1];
outValue += 0.25f * d_r[idX+1 + (idY+1)*n1 + n1+1];
outValue += 0.50f * d_r[idX + (idY-1)*n1 + n1+1];
outValue += 0.50f * d_r[idX + (idY+1)*n1 + n1+1];
这是您询问的 "rows and columns" 方法。这个想法是同时处理输出图像的非重叠区域。
我使用的 16x16 像素示例可以扩展到任何图像尺寸。每个黑框代表一个像素。不同深浅的蓝色仅用于区分要计算的工作项区域。
如果将图像分成两个像素宽的列，则可以在单独的线程中处理每一列。这可以防止列之间的全局写入冲突。对奇数列重复 - 即 dx=1 以计算与列重叠的 2x2 分组，如第二个图像所示。这种方法只是防止了水平方向的写冲突，但是还是需要考虑到行数。
为了避免列内的写入冲突，您需要进一步将列分成行，如上图所示。让我们让一个工作项计算单个 2x2 输出以使内核尽可能简单。行也是分两个阶段计算的——dy=0 和 dy=1。总的来说，这是一个四步算法，其中每一步都可以称为 'embarrassingly parallel'。这些步骤可以通过它们相对于原点的位置偏移来引用：(dx,dy) = A(0,0)、B(1,0)、C(1,0) 和 D(1,1)。请注意，这些组不一定使用相同数量的工作项，具体取决于最终 row/column 是否足以进行计算。
显示内核调用期间两个示例工作项负责的插图。 "rows and columns" 基本上变成了“4 个棋盘”。下面是内核，它将计算 2x2 像素组的平均值并添加全局缓冲区。只要A、B、C、D组各自独立执行，就不需要检查全局写冲突。
_kernel void soSmoothing(__global const float* restrict inData, __global float* restrict outData, int width, int height, int dx, int dy)
{
//(idX,idY) represents the position to compute and write to outData
int idX = get_global_id(0);
int idY = get_global_id(1);
//(pixelX, pixelY) is the top-left corner of the square to compute
int pixelX = idX *2 +dx;
int pixelY = idY *2 +dy;
if(pixelX > (width-2)) return;
if(pixelY > (height-2)) return;
//do the math and add to the corresponding addresses in outData
int topLeftAddress = pixelX + pixelY * width;
float avg = 0.25f * (inData[topLeftAddress] + inData[topLeftAddress +1] + inData[topLeftAddress + width] + inData[topLeftAddress + width +1]);
outData[topLeftAddress] = outData[topLeftAddress] + avg;
outData[topLeftAddress +1] = outData[topLeftAddress +1] + avg;
outData[topLeftAddress +width] = outData[topLeftAddress +width] + avg;
outData[topLeftAddress +width +1] = outData[topLeftAddress +width +1] + avg;
}
主机代码步骤：
1) 将命令队列设置为阻塞
2) 设置缓冲区：1x 输入，1x 输出
3) 将输入缓冲区复制到设备，在设备上将输出初始化为 0
4a) 使用 dx=0, dy=0
使内核入队4b) 使用 dx=1, dy=0
使内核入队4c) 使用 dx=0, dy=1
使内核入队4d) 使用 dx=1, dy=1
使内核入队5) 将输出缓冲区读回主机
请注意，在所有内核执行完成之前，不需要将输出复制到主机。
优点：
- 所有 4 个内核的相同输入 运行 意味着您可以非常轻松地 运行 在最多 4 个独立的设备上执行此操作，并在主机上汇总它们的输出。也许尝试 运行 在 cpu 上安装 1 个内核，在 gpu 上安装 3 个内核。您还可以分解内核以处理图像的较小区域，并根据需要让尽可能多的 opencl 设备参与。传输的开销最终会成为主要瓶颈，因此这可能不值得付出努力。
- 内核简单，易于调试。
- 内核是 'chained'，因此它们会添加到之前的输出中。可以是 运行 任意顺序。
- 内核（+主机代码）可以轻松扩展以支持 3x3 或其他单元尺寸。
缺点：
- 输入图像被读取 4x
- 输出像素写入 4x
- 每个工作项一次只计算 4 个像素。计算更大的区域可以更快，但代码复杂度更高。
给你。请让我知道这对您来说效果如何，以及我是否需要进行任何更正。