CUDA 粒子中的最近邻
Nearest Neighbors in CUDA Particles
编辑 2： 请查看 this crosspost 了解 TLDR。
编辑：鉴于粒子被分割成网格单元（比如16^3网格），让运行一个工作是不是更好的主意- 为每个网格单元分组，一个工作组中的工作项数量与每个网格单元的最大粒子数一样多？
在那种情况下，我可以将相邻单元格中的所有粒子加载到本地内存中，并遍历它们计算一些属性。然后我可以将特定值写入当前网格单元格中的每个粒子。
这种方法是否比 运行 对所有粒子和每次迭代（大部分时间相同）邻居的内核有益？
另外，number of particles/number of grid cells的理想比例是多少？
我正在尝试为 OpenCL 重新实现（和修改）CUDA Particles，并使用它来查询每个粒子的最近邻居。我创建了以下结构：
- 缓冲区
P 保存所有粒子的 3D 位置 (float3)
缓冲区 Sp 存储 int2 对粒子 ID 及其空间哈希。 Sp是根据hash排序的。 （散列只是从 3D 到 1D 的简单线性映射——还没有 Z 索引。）
缓冲区L 在缓冲区Sp 中存储int2 对特定空间散列的起始和结束位置。示例：L[12] = (int2)(0, 50)。
L[12].x 是具有空间散列 12.Sp 中）
L[12].y 是 last 粒子的索引（在 Sp 中），空间散列为 12.
现在我有了所有这些缓冲区，我想遍历 P 中的所有粒子，并为每个粒子遍历其最近的邻居。目前我有一个看起来像这样的内核（伪代码）：
__kernel process_particles(float3* P, int2* Sp, int2* L, int* Out) {
size_t gid = get_global_id(0);
float3 curr_particle = P[gid];
int processed_value = 0;
for(int x=-1; x<=1; x++)
for(int y=-1; y<=1; y++)
for(int z=-1; z<=1; z++) {
float3 neigh_position = curr_particle + (float3)(x,y,z)*GRID_CELL_SIDE;
// ugly boundary checking
if ( dot(neigh_position<0, (float3)(1)) +
dot(neigh_position>BOUNDARY, (float3)(1)) != 0)
continue;
int neigh_hash = spatial_hash( neigh_position );
int2 particles_range = L[ neigh_hash ];
for(int p=particles_range.x; p<particles_range.y; p++)
processed_value += heavy_computation( P[ Sp[p].y ] );
}
Out[gid] = processed_value;
}
该代码的问题在于它很慢。我怀疑非线性 GPU 内存访问（特别是最内层 for 循环中的 P[Sp[p].y]）导致了缓慢。
我想做的是使用Z-order curve作为空间散列。这样我在查询邻居时可以只有 1 for 循环遍历连续的内存范围。唯一的问题是我不知道开始和结束 Z-index 值应该是多少。
我想达到的圣杯：
__kernel process_particles(float3* P, int2* Sp, int2* L, int* Out) {
size_t gid = get_global_id(0);
float3 curr_particle = P[gid];
int processed_value = 0;
// How to accomplish this??
// `get_neighbors_range()` returns start and end Z-index values
// representing the start and end near neighbors cells range
int2 nearest_neighboring_cells_range = get_neighbors_range(curr_particle);
int first_particle_id = L[ nearest_neighboring_cells_range.x ].x;
int last_particle_id = L[ nearest_neighboring_cells_range.y ].y;
for(int p=first_particle_id; p<=last_particle_id; p++) {
processed_value += heavy_computation( P[ Sp[p].y ] );
}
Out[gid] = processed_value;
}
您应该仔细研究 Morton Code 算法。 Ericsons 实时碰撞检测很好地解释了这一点。
Ericson - Real time Collision detection
这是另一个很好的解释，包括一些测试：
Morton encoding/decoding through bit interleaving: Implementations
Z-Order 算法仅定义坐标的路径，您可以在其中从 2 或 3D 坐标散列为整数。尽管算法在每次迭代中都变得更深入，但您必须自己设置限制。通常停止索引由哨兵表示。让哨兵停止将告诉您粒子放置在哪个级别。因此，您要定义的最大级别将告诉您每个维度的单元格数。例如，最大级别为 6 时，您有 2^6 = 64。您的系统 (3D) 中将有 64x64x64 个单元格。这也意味着您必须使用基于整数的坐标。如果您使用浮点数，则必须像 coord.x = 64*float_x 等进行转换。
如果您知道您的系统中有多少个细胞，您就可以定义您的限制。您是否尝试使用二进制八叉树？
由于粒子在运动（在那个 CUDA 示例中），您应该尝试并行化粒子数量而不是细胞数量。
如果你想建立最近邻列表，你必须将粒子映射到细胞。这是通过 table 完成的，然后按细胞对粒子进行分类。你仍然应该遍历粒子并访问它的邻居。
关于您的代码：
The problem with that code is that it's slow. I suspect the nonlinear GPU memory access (particulary P[Sp[p].y] in the inner-most for loop) to be causing the slowness.
记住 Donald Knuth。您应该测量瓶颈的位置。您可以使用 NVCC Profiler 并寻找瓶颈。不确定 OpenCL 有什么分析器。
// ugly boundary checking
if ( dot(neigh_position<0, (float3)(1)) +
dot(neigh_position>BOUNDARY, (float3)(1)) != 0)
continue;
我认为你不应该那样分支，当你调用 heavy_computation 时返回零怎么样？不确定，但也许您在这里有某种分支预测。尝试以某种方式删除它。
运行 只有当您没有对粒子数据的写访问权时，在单元上并行是一个好主意，否则您将不得不使用原子。如果你越过粒子范围而不是读取对单元格和邻居的访问，但你并行创建你的总和并且你不会被迫使用某些竞争条件范例。
Also, what is the ideal ratio of number of particles/number of grid cells?
确实取决于您的算法和您所在域内的粒子堆积，但在您的情况下，我会定义与粒子直径相等的像元大小，并且只使用您获得的像元数。
因此，如果您想使用 Z 顺序并实现您的圣杯，请尝试使用整数坐标并对其进行哈希处理。
也尝试使用更多的粒子。您应该考虑大约 65000 个像 CUDA 示例使用的粒子，因为这样并行化最有效； 运行 处理单元被利用（更少的空闲线程）。
编辑 2： 请查看 this crosspost 了解 TLDR。
编辑：鉴于粒子被分割成网格单元（比如16^3网格），让运行一个工作是不是更好的主意- 为每个网格单元分组，一个工作组中的工作项数量与每个网格单元的最大粒子数一样多？
在那种情况下，我可以将相邻单元格中的所有粒子加载到本地内存中，并遍历它们计算一些属性。然后我可以将特定值写入当前网格单元格中的每个粒子。
这种方法是否比 运行 对所有粒子和每次迭代（大部分时间相同）邻居的内核有益？
另外，number of particles/number of grid cells的理想比例是多少？
我正在尝试为 OpenCL 重新实现（和修改）CUDA Particles，并使用它来查询每个粒子的最近邻居。我创建了以下结构：
- 缓冲区
P保存所有粒子的 3D 位置 (float3) 缓冲区
Sp存储int2对粒子 ID 及其空间哈希。Sp是根据hash排序的。 （散列只是从 3D 到 1D 的简单线性映射——还没有 Z 索引。）缓冲区
L在缓冲区Sp中存储int2对特定空间散列的起始和结束位置。示例：L[12] = (int2)(0, 50)。L[12].x是具有空间散列12. 的 第一个 粒子的索引（在 L[12].y是 last 粒子的索引（在Sp中），空间散列为12.
Sp中）
现在我有了所有这些缓冲区，我想遍历 P 中的所有粒子，并为每个粒子遍历其最近的邻居。目前我有一个看起来像这样的内核（伪代码）：
__kernel process_particles(float3* P, int2* Sp, int2* L, int* Out) {
size_t gid = get_global_id(0);
float3 curr_particle = P[gid];
int processed_value = 0;
for(int x=-1; x<=1; x++)
for(int y=-1; y<=1; y++)
for(int z=-1; z<=1; z++) {
float3 neigh_position = curr_particle + (float3)(x,y,z)*GRID_CELL_SIDE;
// ugly boundary checking
if ( dot(neigh_position<0, (float3)(1)) +
dot(neigh_position>BOUNDARY, (float3)(1)) != 0)
continue;
int neigh_hash = spatial_hash( neigh_position );
int2 particles_range = L[ neigh_hash ];
for(int p=particles_range.x; p<particles_range.y; p++)
processed_value += heavy_computation( P[ Sp[p].y ] );
}
Out[gid] = processed_value;
}
该代码的问题在于它很慢。我怀疑非线性 GPU 内存访问（特别是最内层 for 循环中的 P[Sp[p].y]）导致了缓慢。
我想做的是使用Z-order curve作为空间散列。这样我在查询邻居时可以只有 1 for 循环遍历连续的内存范围。唯一的问题是我不知道开始和结束 Z-index 值应该是多少。
我想达到的圣杯：
__kernel process_particles(float3* P, int2* Sp, int2* L, int* Out) {
size_t gid = get_global_id(0);
float3 curr_particle = P[gid];
int processed_value = 0;
// How to accomplish this??
// `get_neighbors_range()` returns start and end Z-index values
// representing the start and end near neighbors cells range
int2 nearest_neighboring_cells_range = get_neighbors_range(curr_particle);
int first_particle_id = L[ nearest_neighboring_cells_range.x ].x;
int last_particle_id = L[ nearest_neighboring_cells_range.y ].y;
for(int p=first_particle_id; p<=last_particle_id; p++) {
processed_value += heavy_computation( P[ Sp[p].y ] );
}
Out[gid] = processed_value;
}
您应该仔细研究 Morton Code 算法。 Ericsons 实时碰撞检测很好地解释了这一点。
Ericson - Real time Collision detection
这是另一个很好的解释，包括一些测试：
Morton encoding/decoding through bit interleaving: Implementations
Z-Order 算法仅定义坐标的路径，您可以在其中从 2 或 3D 坐标散列为整数。尽管算法在每次迭代中都变得更深入，但您必须自己设置限制。通常停止索引由哨兵表示。让哨兵停止将告诉您粒子放置在哪个级别。因此，您要定义的最大级别将告诉您每个维度的单元格数。例如，最大级别为 6 时，您有 2^6 = 64。您的系统 (3D) 中将有 64x64x64 个单元格。这也意味着您必须使用基于整数的坐标。如果您使用浮点数，则必须像 coord.x = 64*float_x 等进行转换。
如果您知道您的系统中有多少个细胞，您就可以定义您的限制。您是否尝试使用二进制八叉树？
由于粒子在运动（在那个 CUDA 示例中），您应该尝试并行化粒子数量而不是细胞数量。
如果你想建立最近邻列表，你必须将粒子映射到细胞。这是通过 table 完成的，然后按细胞对粒子进行分类。你仍然应该遍历粒子并访问它的邻居。
关于您的代码：
The problem with that code is that it's slow. I suspect the nonlinear GPU memory access (particulary P[Sp[p].y] in the inner-most for loop) to be causing the slowness.
记住 Donald Knuth。您应该测量瓶颈的位置。您可以使用 NVCC Profiler 并寻找瓶颈。不确定 OpenCL 有什么分析器。
// ugly boundary checking if ( dot(neigh_position<0, (float3)(1)) + dot(neigh_position>BOUNDARY, (float3)(1)) != 0) continue;
我认为你不应该那样分支，当你调用 heavy_computation 时返回零怎么样？不确定，但也许您在这里有某种分支预测。尝试以某种方式删除它。
运行 只有当您没有对粒子数据的写访问权时，在单元上并行是一个好主意，否则您将不得不使用原子。如果你越过粒子范围而不是读取对单元格和邻居的访问，但你并行创建你的总和并且你不会被迫使用某些竞争条件范例。
Also, what is the ideal ratio of number of particles/number of grid cells?
确实取决于您的算法和您所在域内的粒子堆积，但在您的情况下，我会定义与粒子直径相等的像元大小，并且只使用您获得的像元数。
因此，如果您想使用 Z 顺序并实现您的圣杯，请尝试使用整数坐标并对其进行哈希处理。
也尝试使用更多的粒子。您应该考虑大约 65000 个像 CUDA 示例使用的粒子，因为这样并行化最有效； 运行 处理单元被利用（更少的空闲线程）。