CUDA 缩减：变形展开（学校）
CUDA Reduction: Warp Unrolling (School)
我目前正在开展一个项目，在该项目中展开最后一次减少。我已经完成了上面的代码；但是，有些修改是通过猜测完成的，我想解释一下原因。我写的代码只有函数kernel4
// in is input array, out is where to store result, n is number of elements from in
// T is a float (32bit)
__global__ void kernel4(T *in, T *out, unsigned int n)
这是一个归约算法，其余代码已经提供。
代码：
#include <stdlib.h>
#include <stdio.h>
#include "timer.h"
#include "cuda_utils.h"
typedef float T;
#define N_ (8 * 1024 * 1024)
#define MAX_THREADS 256
#define MAX_BLOCKS 64
#define MIN(x,y) ((x < y) ? x : y)
#define tid threadIdx.x
#define bid blockIdx.x
#define bdim blockDim.x
#define warp_size 32
unsigned int nextPow2( unsigned int x ) {
--x;
x |= x >> 1;
x |= x >> 2;
x |= x >> 4;
x |= x >> 8;
x |= x >> 16;
return ++x;
}
void getNumBlocksAndThreads(int whichKernel, int n, int maxBlocks, int maxThreads, int &blocks, int &threads)
{
if (whichKernel < 3) {
threads = (n < maxThreads) ? nextPow2(n) : maxThreads;
blocks = (n + threads - 1) / threads;
} else {
threads = (n < maxThreads*2) ? nextPow2((n + 1)/ 2) : maxThreads;
blocks = (n + (threads * 2 - 1)) / (threads * 2);
}
if (whichKernel == 5)
blocks = MIN(maxBlocks, blocks);
}
T reduce_cpu(T *data, int n) {
T sum = data[0];
T c = (T) 0.0;
for (int i = 1; i < n; i++)
{
T y = data[i] - c;
T t = sum + y;
c = (t - sum) - y;
sum = t;
}
return sum;
}
__global__ void
kernel4(T *in, T *out, unsigned int n)
{
__shared__ volatile T d[MAX_THREADS];
unsigned int i = bid * bdim + tid;
n >>= 1;
d[tid] = (i < n) ? in[i] + in[i+n] : 0;
__syncthreads ();
for(unsigned int s = bdim >> 1; s > warp_size; s >>= 1) {
if(tid < s)
d[tid] += d[tid + s];
__syncthreads ();
}
if (tid < warp_size) {
if (n > 64) d[tid] += d[tid + 32];
if (n > 32) d[tid] += d[tid + 16];
d[tid] += d[tid + 8];
d[tid] += d[tid + 4];
d[tid] += d[tid + 2];
d[tid] += d[tid + 1];
}
if(tid == 0)
out[bid] = d[0];
}
int main(int argc, char** argv)
{
T *h_idata, h_odata, h_cpu;
T *d_idata, *d_odata;
struct stopwatch_t* timer = NULL;
long double t_kernel_4, t_cpu;
int whichKernel = 4, threads, blocks, N, i;
if(argc > 1) {
N = atoi (argv[1]);
printf("N: %d\n", N);
} else {
N = N_;
printf("N: %d\n", N);
}
getNumBlocksAndThreads (whichKernel, N, MAX_BLOCKS, MAX_THREADS, blocks, threads);
stopwatch_init ();
timer = stopwatch_create ();
h_idata = (T*) malloc (N * sizeof (T));
CUDA_CHECK_ERROR (cudaMalloc (&d_idata, N * sizeof (T)));
CUDA_CHECK_ERROR (cudaMalloc (&d_odata, blocks * sizeof (T)));
srand48(time(NULL));
for(i = 0; i < N; i++)
h_idata[i] = drand48() / 100000;
CUDA_CHECK_ERROR (cudaMemcpy (d_idata, h_idata, N * sizeof (T), cudaMemcpyHostToDevice));
dim3 gb(blocks, 1, 1);
dim3 tb(threads, 1, 1);
kernel4 <<<gb, tb>>> (d_idata, d_odata, N);
cudaThreadSynchronize ();
stopwatch_start (timer);
kernel4 <<<gb, tb>>> (d_idata, d_odata, N);
int s = blocks;
while(s > 1) {
threads = 0;
blocks = 0;
getNumBlocksAndThreads (whichKernel, s, MAX_BLOCKS, MAX_THREADS, blocks, threads);
dim3 gb(blocks, 1, 1);
dim3 tb(threads, 1, 1);
kernel4 <<<gb, tb>>> (d_odata, d_odata, s);
s = (s + threads * 2 - 1) / (threads * 2);
}
cudaThreadSynchronize ();
t_kernel_4 = stopwatch_stop (timer);
fprintf (stdout, "Time to execute unrolled GPU reduction kernel: %Lg secs\n", t_kernel_4);
double bw = (N * sizeof(T)) / (t_kernel_4 * 1e9); // total bits / time
fprintf (stdout, "Effective bandwidth: %.2lf GB/s\n", bw);
CUDA_CHECK_ERROR (cudaMemcpy (&h_odata, d_odata, sizeof (T), cudaMemcpyDeviceToHost));
stopwatch_start (timer);
h_cpu = reduce_cpu (h_idata, N);
t_cpu = stopwatch_stop (timer);
fprintf (stdout, "Time to execute naive CPU reduction: %Lg secs\n", t_cpu);
if(abs (h_odata - h_cpu) > 1e-5)
fprintf(stderr, "FAILURE: GPU: %f CPU: %f\n", h_odata, h_cpu);
else
printf("SUCCESS: GPU: %f CPU: %f\n", h_odata, h_cpu);
return 0;
}
我的第一个问题是：声明时
__shared__ volatile T d[MAX_THREADS];
我想验证一下我对volatile的理解。 Volatile 防止编译器错误地优化我的代码，并承诺 load/stores 是通过缓存完成的，而不仅仅是寄存器（如果错误请纠正我）。对于归约，如果部分归约和仍然存储在寄存器中，为什么会出现问题？
我的第二个问题是：在进行实际的 warp reduction 时
if (tid < warp_size) { // Final log2(32) = 5 strides
if (n > 64) d[tid] += d[tid + 32];
if (n > 32) d[tid] += d[tid + 16];
d[tid] += d[tid + 8];
d[tid] += d[tid + 4];
d[tid] += d[tid + 2];
d[tid] += d[tid + 1];
}
如果没有 (n > 64) 和 (n > 32) 条件，归约和将产生不正确的结果。我得到的结果是：
FAILURE: GPU: 41.966557 CPU: 41.946209
经过 5 次试验，GPU 缩减始终产生 0.0204 的误差。我谨慎地认为这是一个浮点运算错误。
老实说，我的老师的助手建议进行此更改以添加 (n > 64) 和 (n > 32) 条件，但没有解释为什么它会修复代码。
既然我的试验中的 n 都超过了 64，为什么这个条件会改变结果。我很难追溯问题，因为我无法像在 CPU.
中那样使用打印功能
在我们解决您的两个问题之前，让我们从一些序言评论开始：
- 我鼓励您阅读 NVIDIA 的 canonical reduction tutorial
- 像这样写的缩减做了几个假设，其中之一是块大小是 power-of-2（对于 "correctness"）。
- 您的代码在最后的缩减阶段使用 warp-synchronous 编程。你似乎知道你在做什么，所以我不会提供详细的描述，但它肯定与理解这里相关。您可以 google 它并在需要时获取说明。它与下面的讨论相关，但我不会在每种情况下都指出它的相关性。
好的，现在你的问题：
I would like to verify my understanding of volatile. Volatile prevents compilers from incorrectly optimizing my code and promises that load/stores are completed through the cache and not just registers (please correct me if wrong). For reduction, if partial reduction sums are still stored in registers, why is this a problem?
关于volatile的定义，请参考the CUDA programming guide。我已经看到摘要描述将此称为阻止寄存器优化或阻止加载和存储的重新排序。我更喜欢前者并将其用作工作定义。
基本思想是 volatile 强制对该变量的任何引用（读或写）实际上进入内存子系统。我的意思是它将执行读取或写入，并且不会尝试使用先前加载到寄存器中的值。如果没有此限定符，编译器可以（例如）从实际内存位置自由加载一个值，然后在它认为合适的时间内将该值（及其任何更新）保存在寄存器中。编译器会着眼于性能来做到这一点。 （顺便说一句，请注意您在这里使用了 "cache" 这个词。我会在这里避免这种用法。共享内存在它和处理器 load/store 机制之间没有插入缓存。）
在这种类型的 warp-synchronous 编码中没有 volatile，如果我们允许编译器 "optimize"（即维护）中间值，我们将 运行 变成一个问题寄存器。这主要是由于inter-thread 通信。为了清楚地了解原因，让我们看看最终减少的最后 2 个步骤：
d[tid] += d[tid + 2];
d[tid] += d[tid + 1];
让我们只考虑 tid 值为 0-1 的线程。在second-last步骤中，线程0将获取d[2]值并将其添加到d[0]值，而线程1将获取d[3]值并将其添加到d[1] 值。此时，如果我们不使用 volatile，编译器就没有义务将线程 1 累积的 d[1] 值写回共享内存。允许将其保存在寄存器中。所以在共享内存中看到的 d[1] 值不是 "up-to-date".
现在让我们进入最后一步。在这一步中，线程0从共享内存中读取d[1]值，并将其添加到d[0]值中。但是如果没有volatile，我们在前面的步骤中看到d[1]的共享内存内容不再准确。 OTOH，如果我们使用 volatile，那么上一步中对共享内存的写入将实际发生，并且在最后一步中，线程 0 将在读取 d[1] 时获取正确的值。 CUDA 线程是一个独立的模型。我的意思是，一个线程不能直接访问属于另一个线程的寄存器中包含的值。因此 inter-thread warp 级别的通信通常通过共享内存或通过 warp-shuffle 操作来完成。
__syncthreads() 有类似的行为：它强制所有像这样的 register-optimized 值被写出到内存，以便它们 "visible" 到块中的其他线程。因此，更复杂的优化是仅在归约从基于 loop-driven __syncthreads() 的归约切换到最终的 warp-synchronous 归约时切换到 volatile 限定指针。您可以在我在此答案开头链接的教程幻灯片中看到一个示例。
另外，warp-synchronous 这种编程在 CUDA 9 中（更正式地）已弃用。相反，您应该使用 cooperative groups.
The reduction sum will yield incorrect results without (n > 64) and (n > 32) conditions.
主要使用这些条件是因为代码被设计为 "correct" 用于具有 power-of-2 大小的任何块配置。如果我们假设块大小（每个块的线程数）是 2 的幂，并且大于 64，则它必须是 128 或更大，例如。您的 n 变量以块大小开始，但随后乘以 2：
n >>= 1;
因此，如果我们要保证这行代码的正确性：
d[tid] += d[tid + 32];
那么我们应该只在线程块大小为 64（至少）时应用该操作，这就像说 n 大于 64：
if (n > 64) d[tid] += d[tid + 32];
关于这个问题，声称如果包含或不包含 if (n > 64)，发布的代码的行为会有所不同。这样做的原因是发布的代码包含一个循环，该循环在减少过程中重新计算线程数和块数：
int s = blocks;
while(s > 1) {
threads = 0;
blocks = 0;
getNumBlocksAndThreads (whichKernel, s, MAX_BLOCKS, MAX_THREADS, blocks, threads);
此循环最终导致块大小小于 128，这意味着省略 if 条件会导致中断。 （在这个循环中简单地打印出 threads 变量）
对此：
I am having difficulty tracing back the problem because I cannot use print functions like I would in a CPU.
我不确定那里有什么问题。 printf 应该在内核代码中工作。
共享 变量不能将初始化作为其声明的一部分 according to this answer。
因此，如果 n < 64，我们将一些随机共享内存数组数据添加到总和中，这种情况下会出错。
我目前正在开展一个项目，在该项目中展开最后一次减少。我已经完成了上面的代码；但是，有些修改是通过猜测完成的，我想解释一下原因。我写的代码只有函数kernel4
// in is input array, out is where to store result, n is number of elements from in
// T is a float (32bit)
__global__ void kernel4(T *in, T *out, unsigned int n)
这是一个归约算法，其余代码已经提供。
代码：
#include <stdlib.h>
#include <stdio.h>
#include "timer.h"
#include "cuda_utils.h"
typedef float T;
#define N_ (8 * 1024 * 1024)
#define MAX_THREADS 256
#define MAX_BLOCKS 64
#define MIN(x,y) ((x < y) ? x : y)
#define tid threadIdx.x
#define bid blockIdx.x
#define bdim blockDim.x
#define warp_size 32
unsigned int nextPow2( unsigned int x ) {
--x;
x |= x >> 1;
x |= x >> 2;
x |= x >> 4;
x |= x >> 8;
x |= x >> 16;
return ++x;
}
void getNumBlocksAndThreads(int whichKernel, int n, int maxBlocks, int maxThreads, int &blocks, int &threads)
{
if (whichKernel < 3) {
threads = (n < maxThreads) ? nextPow2(n) : maxThreads;
blocks = (n + threads - 1) / threads;
} else {
threads = (n < maxThreads*2) ? nextPow2((n + 1)/ 2) : maxThreads;
blocks = (n + (threads * 2 - 1)) / (threads * 2);
}
if (whichKernel == 5)
blocks = MIN(maxBlocks, blocks);
}
T reduce_cpu(T *data, int n) {
T sum = data[0];
T c = (T) 0.0;
for (int i = 1; i < n; i++)
{
T y = data[i] - c;
T t = sum + y;
c = (t - sum) - y;
sum = t;
}
return sum;
}
__global__ void
kernel4(T *in, T *out, unsigned int n)
{
__shared__ volatile T d[MAX_THREADS];
unsigned int i = bid * bdim + tid;
n >>= 1;
d[tid] = (i < n) ? in[i] + in[i+n] : 0;
__syncthreads ();
for(unsigned int s = bdim >> 1; s > warp_size; s >>= 1) {
if(tid < s)
d[tid] += d[tid + s];
__syncthreads ();
}
if (tid < warp_size) {
if (n > 64) d[tid] += d[tid + 32];
if (n > 32) d[tid] += d[tid + 16];
d[tid] += d[tid + 8];
d[tid] += d[tid + 4];
d[tid] += d[tid + 2];
d[tid] += d[tid + 1];
}
if(tid == 0)
out[bid] = d[0];
}
int main(int argc, char** argv)
{
T *h_idata, h_odata, h_cpu;
T *d_idata, *d_odata;
struct stopwatch_t* timer = NULL;
long double t_kernel_4, t_cpu;
int whichKernel = 4, threads, blocks, N, i;
if(argc > 1) {
N = atoi (argv[1]);
printf("N: %d\n", N);
} else {
N = N_;
printf("N: %d\n", N);
}
getNumBlocksAndThreads (whichKernel, N, MAX_BLOCKS, MAX_THREADS, blocks, threads);
stopwatch_init ();
timer = stopwatch_create ();
h_idata = (T*) malloc (N * sizeof (T));
CUDA_CHECK_ERROR (cudaMalloc (&d_idata, N * sizeof (T)));
CUDA_CHECK_ERROR (cudaMalloc (&d_odata, blocks * sizeof (T)));
srand48(time(NULL));
for(i = 0; i < N; i++)
h_idata[i] = drand48() / 100000;
CUDA_CHECK_ERROR (cudaMemcpy (d_idata, h_idata, N * sizeof (T), cudaMemcpyHostToDevice));
dim3 gb(blocks, 1, 1);
dim3 tb(threads, 1, 1);
kernel4 <<<gb, tb>>> (d_idata, d_odata, N);
cudaThreadSynchronize ();
stopwatch_start (timer);
kernel4 <<<gb, tb>>> (d_idata, d_odata, N);
int s = blocks;
while(s > 1) {
threads = 0;
blocks = 0;
getNumBlocksAndThreads (whichKernel, s, MAX_BLOCKS, MAX_THREADS, blocks, threads);
dim3 gb(blocks, 1, 1);
dim3 tb(threads, 1, 1);
kernel4 <<<gb, tb>>> (d_odata, d_odata, s);
s = (s + threads * 2 - 1) / (threads * 2);
}
cudaThreadSynchronize ();
t_kernel_4 = stopwatch_stop (timer);
fprintf (stdout, "Time to execute unrolled GPU reduction kernel: %Lg secs\n", t_kernel_4);
double bw = (N * sizeof(T)) / (t_kernel_4 * 1e9); // total bits / time
fprintf (stdout, "Effective bandwidth: %.2lf GB/s\n", bw);
CUDA_CHECK_ERROR (cudaMemcpy (&h_odata, d_odata, sizeof (T), cudaMemcpyDeviceToHost));
stopwatch_start (timer);
h_cpu = reduce_cpu (h_idata, N);
t_cpu = stopwatch_stop (timer);
fprintf (stdout, "Time to execute naive CPU reduction: %Lg secs\n", t_cpu);
if(abs (h_odata - h_cpu) > 1e-5)
fprintf(stderr, "FAILURE: GPU: %f CPU: %f\n", h_odata, h_cpu);
else
printf("SUCCESS: GPU: %f CPU: %f\n", h_odata, h_cpu);
return 0;
}
我的第一个问题是：声明时
__shared__ volatile T d[MAX_THREADS];
我想验证一下我对volatile的理解。 Volatile 防止编译器错误地优化我的代码，并承诺 load/stores 是通过缓存完成的，而不仅仅是寄存器（如果错误请纠正我）。对于归约，如果部分归约和仍然存储在寄存器中，为什么会出现问题？
我的第二个问题是：在进行实际的 warp reduction 时
if (tid < warp_size) { // Final log2(32) = 5 strides
if (n > 64) d[tid] += d[tid + 32];
if (n > 32) d[tid] += d[tid + 16];
d[tid] += d[tid + 8];
d[tid] += d[tid + 4];
d[tid] += d[tid + 2];
d[tid] += d[tid + 1];
}
如果没有 (n > 64) 和 (n > 32) 条件，归约和将产生不正确的结果。我得到的结果是：
FAILURE: GPU: 41.966557 CPU: 41.946209
经过 5 次试验，GPU 缩减始终产生 0.0204 的误差。我谨慎地认为这是一个浮点运算错误。
老实说，我的老师的助手建议进行此更改以添加 (n > 64) 和 (n > 32) 条件，但没有解释为什么它会修复代码。
既然我的试验中的 n 都超过了 64，为什么这个条件会改变结果。我很难追溯问题，因为我无法像在 CPU.
中那样使用打印功能在我们解决您的两个问题之前，让我们从一些序言评论开始：
- 我鼓励您阅读 NVIDIA 的 canonical reduction tutorial
- 像这样写的缩减做了几个假设，其中之一是块大小是 power-of-2（对于 "correctness"）。
- 您的代码在最后的缩减阶段使用 warp-synchronous 编程。你似乎知道你在做什么，所以我不会提供详细的描述，但它肯定与理解这里相关。您可以 google 它并在需要时获取说明。它与下面的讨论相关，但我不会在每种情况下都指出它的相关性。
好的，现在你的问题：
I would like to verify my understanding of volatile. Volatile prevents compilers from incorrectly optimizing my code and promises that load/stores are completed through the cache and not just registers (please correct me if wrong). For reduction, if partial reduction sums are still stored in registers, why is this a problem?
关于volatile的定义，请参考the CUDA programming guide。我已经看到摘要描述将此称为阻止寄存器优化或阻止加载和存储的重新排序。我更喜欢前者并将其用作工作定义。
基本思想是 volatile 强制对该变量的任何引用（读或写）实际上进入内存子系统。我的意思是它将执行读取或写入，并且不会尝试使用先前加载到寄存器中的值。如果没有此限定符，编译器可以（例如）从实际内存位置自由加载一个值，然后在它认为合适的时间内将该值（及其任何更新）保存在寄存器中。编译器会着眼于性能来做到这一点。 （顺便说一句，请注意您在这里使用了 "cache" 这个词。我会在这里避免这种用法。共享内存在它和处理器 load/store 机制之间没有插入缓存。）
在这种类型的 warp-synchronous 编码中没有 volatile，如果我们允许编译器 "optimize"（即维护）中间值，我们将 运行 变成一个问题寄存器。这主要是由于inter-thread 通信。为了清楚地了解原因，让我们看看最终减少的最后 2 个步骤：
d[tid] += d[tid + 2];
d[tid] += d[tid + 1];
让我们只考虑 tid 值为 0-1 的线程。在second-last步骤中，线程0将获取d[2]值并将其添加到d[0]值，而线程1将获取d[3]值并将其添加到d[1] 值。此时，如果我们不使用 volatile，编译器就没有义务将线程 1 累积的 d[1] 值写回共享内存。允许将其保存在寄存器中。所以在共享内存中看到的 d[1] 值不是 "up-to-date".
现在让我们进入最后一步。在这一步中，线程0从共享内存中读取d[1]值，并将其添加到d[0]值中。但是如果没有volatile，我们在前面的步骤中看到d[1]的共享内存内容不再准确。 OTOH，如果我们使用 volatile，那么上一步中对共享内存的写入将实际发生，并且在最后一步中，线程 0 将在读取 d[1] 时获取正确的值。 CUDA 线程是一个独立的模型。我的意思是，一个线程不能直接访问属于另一个线程的寄存器中包含的值。因此 inter-thread warp 级别的通信通常通过共享内存或通过 warp-shuffle 操作来完成。
__syncthreads() 有类似的行为：它强制所有像这样的 register-optimized 值被写出到内存，以便它们 "visible" 到块中的其他线程。因此，更复杂的优化是仅在归约从基于 loop-driven __syncthreads() 的归约切换到最终的 warp-synchronous 归约时切换到 volatile 限定指针。您可以在我在此答案开头链接的教程幻灯片中看到一个示例。
另外，warp-synchronous 这种编程在 CUDA 9 中（更正式地）已弃用。相反，您应该使用 cooperative groups.
The reduction sum will yield incorrect results without (n > 64) and (n > 32) conditions.
主要使用这些条件是因为代码被设计为 "correct" 用于具有 power-of-2 大小的任何块配置。如果我们假设块大小（每个块的线程数）是 2 的幂，并且大于 64，则它必须是 128 或更大，例如。您的 n 变量以块大小开始，但随后乘以 2：
n >>= 1;
因此，如果我们要保证这行代码的正确性：
d[tid] += d[tid + 32];
那么我们应该只在线程块大小为 64（至少）时应用该操作，这就像说 n 大于 64：
if (n > 64) d[tid] += d[tid + 32];
关于这个问题，声称如果包含或不包含 if (n > 64)，发布的代码的行为会有所不同。这样做的原因是发布的代码包含一个循环，该循环在减少过程中重新计算线程数和块数：
int s = blocks;
while(s > 1) {
threads = 0;
blocks = 0;
getNumBlocksAndThreads (whichKernel, s, MAX_BLOCKS, MAX_THREADS, blocks, threads);
此循环最终导致块大小小于 128，这意味着省略 if 条件会导致中断。 （在这个循环中简单地打印出 threads 变量）
对此：
I am having difficulty tracing back the problem because I cannot use print functions like I would in a CPU.
我不确定那里有什么问题。 printf 应该在内核代码中工作。
共享 变量不能将初始化作为其声明的一部分 according to this answer。 因此，如果 n < 64，我们将一些随机共享内存数组数据添加到总和中，这种情况下会出错。