你在删除的答案中得出的主要结论是正确的：你发布的内核没有意识到在内核执行结束时，你已经完成了整体减少的大部分工作，但结果并不完全。每个块的结果必须被组合起来（以某种方式）。正如评论中指出的那样，你的代码还有一些其他问题。让我们看看修改后的版本：
__device__ float atomicMaxf(float* address, float val)
{
int *address_as_int =(int*)address;
int old = *address_as_int, assumed;
while (val > __int_as_float(old)) {
assumed = old;
old = atomicCAS(address_as_int, assumed,
__float_as_int(val));
}
return __int_as_float(old);
}
__global__ void max_reduce(const float* const d_array, float* d_max,
const size_t elements)
{
extern __shared__ float shared[];
int tid = threadIdx.x;
int gid = (blockDim.x * blockIdx.x) + tid;
shared[tid] = -FLOAT_MAX;
if (gid < elements)
shared[tid] = d_array[gid];
__syncthreads();
for (unsigned int s=blockDim.x/2; s>0; s>>=1)
{
if (tid < s && gid < elements)
shared[tid] = max(shared[tid], shared[tid + s]);
__syncthreads();
}
if (tid == 0)
d_max[blockIdx.x] = shared[tid];
if (tid == 0)
atomicMaxf(d_max, shared[0]);
}
- 正如Pavan所述，您需要初始化共享内存数组。 如果
gridDim.x * blockDim.x大于elements，则最后启动的块可能不是“完整”块。
- 请注意，在此行中，尽管我们正在检查操作线程（
gid）是否小于elements，但当我们将添加到gid用于索引进入共享内存，我们仍然可以在最后一个块中超出合法值复制到共享内存之外进行索引。 因此，我们需要注意1中指示的共享内存初始化。
- 正如您已经发现的那样，您的最后一行是不正确的。 每个块都会产生自己的结果，我们必须以某种方式将它们组合起来。 如果启动的块数较小（稍后会详细说明），您可以考虑使用原子操作之一。 通常我们不建议人们使用原子操作，因为它们在执行时间方面是“昂贵”的。 但是，我们面临的另一个选择是将块结果保存在全局内存中，完成内核，然后可能启动另一个内核来合并各个块结果。 如果我最初启动了大量块（比如超过1024个），如果我按照这种方法进行操作，我可能最终会启动两个附加内核。 因此考虑原子操作。 如上所示，对于浮点数没有本地
atomicMax函数，但是正如文档中所述，您可以使用atomicCAS生成任何任意的原子函数，并且我已经在atomicMaxf中提供了一个示例，该示例为float提供原子最大值。
但是运行1024个或更多个原子函数（每个块一个）是否是最好的方法？ 可能不是。
在启动线程块的内核时，我们实际上只需要启动足够的线程块来使机器保持繁忙状态。 作为经验法则，我们希望每个SM上至少有4-8个warp正在操作，稍微多一些可能是一个好主意。 但是从机器利用率的角度来看，最初启动数千个线程块并没有特别的好处。 如果我们选择像每个SM 8个线程块这样的数字，而我们的GPU最多只有14-16个SM，这给我们了一个相对较小的线程块数8 * 14 = 112个。 让我们选择128（8 * 16）作为一个漂亮的圆整数字。 这并没有什么神奇之处，只是足以使GPU保持繁忙状态。 如果我们让这128个线程块做额外的工作来解决整个问题，那么我们就可以利用原子操作而不必支付太高的代价，并避免多次内核启动。 那么这会是什么样子呢？
__device__ float atomicMaxf(float* address, float val)
{
int *address_as_int =(int*)address;
int old = *address_as_int, assumed;
while (val > __int_as_float(old)) {
assumed = old;
old = atomicCAS(address_as_int, assumed,
__float_as_int(val));
}
return __int_as_float(old);
}
__global__ void max_reduce(const float* const d_array, float* d_max,
const size_t elements)
{
extern __shared__ float shared[];
int tid = threadIdx.x;
int gid = (blockDim.x * blockIdx.x) + tid;
shared[tid] = -FLOAT_MAX;
while (gid < elements) {
shared[tid] = max(shared[tid], d_array[gid]);
gid += gridDim.x*blockDim.x;
}
__syncthreads();
gid = (blockDim.x * blockIdx.x) + tid;
for (unsigned int s=blockDim.x/2; s>0; s>>=1)
{
if (tid < s && gid < elements)
shared[tid] = max(shared[tid], shared[tid + s]);
__syncthreads();
}
if (tid == 0)
atomicMaxf(d_max, shared[0]);
}
使用这个修改过的内核，在创建内核启动时，我们不是根据整个数据大小（elements）来决定启动多少个线程块，而是启动一定数量的块（比如说128个，您可以修改此数字以找出运行最快的），并让每个线程块（因此整个网格）在共享内存中循环，对每个元素进行部分最大操作。然后，在标有注释1的行中，我们必须将gid变量重新设置为其初始值。这实际上是不必要的，如果我们保证网格的大小（gridDim.x*blockDim.x）小于elements，则可以进一步简化块约简循环代码，在内核启动时很容易做到这一点。
请注意，使用这种原子方法时，有必要将结果（在这种情况下是*d_max）初始化为适当的值，如-FLOAT_MAX。
再次说明，通常我们不建议使用原子操作，但在这种情况下，如果我们仔细管理它，并且它允许我们节省额外的内核启动开销，那么考虑使用它是值得的。
要了解如何进行快速并行约简的忍者级别分析，请查看Mark Harris的优秀白皮书，该白皮书可与相关CUDA示例一起使用。
