CUDA并行前缀和错误
CUDA Paralell prefix sum error
我正在尝试按照《大规模并行处理器编程》第 3 版第 8 章中的描述实现三阶段并行扫描(有任何代码行,但只有指令)。
该算法只允许使用 1 个块,块中线程数最大,并且受共享内存大小的限制,因为所有元素都必须适合共享内存
经过一些调试,我在第 3 阶段求和时遇到了问题(见下面的代码),当我使用大量元素,例如 8192,并且超过 1 个线程时。
算法的图形概念如下:
下面是内核代码:
__global__
void efficient_Kogge_Stone_scan_kernel(float *X, float *Y, int InputSize) {
__shared__ float XY[SECTION_SIZE];
__shared__ float AUS[BLOCK_DIM];
//int i = blockIdx.x * blockDim.x + threadIdx.x;
// Keep mind: Partition the input into blockDim.x subsections: i.e. for 8 threads --> 8 subsections
// collaborative load in a coalesced manner
for (int j = 0; j < SECTION_SIZE; j += blockDim.x) {
XY[threadIdx.x + j] = X[threadIdx.x + j];
}
__syncthreads();
// PHASE 1: scan inner own subsection
// At the end of this phase the last element of each subsection contains the sum of all alements in own subsection
for (int j = 1; j < SUBSECTION_SIZE; j++) {
XY[threadIdx.x * (SUBSECTION_SIZE)+j] += XY[threadIdx.x * (SUBSECTION_SIZE)+j - 1];
}
__syncthreads();
// PHASE 2: perform iterative kogge_stone_scan of the last elements of each subsections of XY loaded first in AUS
AUS[threadIdx.x] = XY[threadIdx.x * (SUBSECTION_SIZE)+(SUBSECTION_SIZE)-1];
float in;
for (unsigned int stride = 1; stride < blockDim.x; stride *= 2) {
__syncthreads();
if (threadIdx.x >= stride) {
in = AUS[threadIdx.x - stride];
}
__syncthreads();
if (threadIdx.x >= stride) {
AUS[threadIdx.x] += in;
}
}
__syncthreads();
// PHASE 3: each thread adds to its elements the new value of the last element of its predecessor's section
if (threadIdx.x > 0) {
for (unsigned int stride = 0; stride < (SUBSECTION_SIZE); stride++) {
XY[threadIdx.x * (SUBSECTION_SIZE)+stride] += AUS[threadIdx.x - 1]; // <--
}
}
__syncthreads();
// store the result into output vector
for (int j = 0; j < SECTION_SIZE; j += blockDim.x) {
Y[threadIdx.x + j] = XY[threadIdx.x + j];
}
}
如果我在块中使用一个线程和 8192 个元素,它会完美地工作,但如果我使用多个线程,则 XY[5793](或 X[5793] 完成并存储时的结果是错误的结果)。
它可以使用 4096 个元素和一个或多个线程(最多 1024 个线程)。
如果我使用 int 而不是 float 数字,它甚至可以使用一个或多个线程的 8192 个元素。
我也尝试在 MATLAB 中进行验证,这些是输出比较:
- X[5973] = 16788115 ---- MATLAB
- X[5973] = 16788114 ---- CPU
- X[5973] = 16788116 ---- GPU
正如我们所见,CPU 结果也与 MATLAB 不同,所以根据这些结果我认为问题出在浮点加法上,但我告诉你我用有序 "x.00" 浮点数(例如 {1.00, 2.00, 3.00, 4.00 ..... 8192.00})。
另一个问题是关于性能,主机代码总是比内核代码快,有这些配置参数和这些输入,正常吗?
如果你需要完整的源代码,你可以找到它here
第一次扫描可能有问题:
XY[threadIdx.x * (SUBSECTION_SIZE)+j] += XY[threadIdx.x * (SUBSECTION_SIZE)+j - 1];
这可能导致对共享内存中元素的读取不连贯。当您读取前一个元素时,不能保证任何其他线程都没有更新该值。
尝试通过将值存储在寄存器中来将此部分分成几部分。示例:
int t = XY[threadIdx.x * (SUBSECTION_SIZE)+j - 1];
__syncthreads();
XY[threadIdx.x * (SUBSECTION_SIZE)+j] += t;
8192 是 2^13
sum(1..8192) 接近8192^2/2: 8192*8193/2,比2^25多一点。
因此您需要 26 位来表示它(请参阅下面的注释)。
单精度 IEEE 754 浮点数只有 24 位有效数,因此,取决于求和的执行方式(以何种顺序),以及最终的舍入方向(通常它是默认舍入到最近的,绑定到偶数),然后结果可能不同。
注意严格来说,精确的和可以用float来表示,不用四舍五入,因为最后12位是0,所以尾数只跨越14位。但部分和不是这样。
我正在尝试按照《大规模并行处理器编程》第 3 版第 8 章中的描述实现三阶段并行扫描(有任何代码行,但只有指令)。 该算法只允许使用 1 个块,块中线程数最大,并且受共享内存大小的限制,因为所有元素都必须适合共享内存
经过一些调试,我在第 3 阶段求和时遇到了问题(见下面的代码),当我使用大量元素,例如 8192,并且超过 1 个线程时。
算法的图形概念如下:
下面是内核代码:
__global__
void efficient_Kogge_Stone_scan_kernel(float *X, float *Y, int InputSize) {
__shared__ float XY[SECTION_SIZE];
__shared__ float AUS[BLOCK_DIM];
//int i = blockIdx.x * blockDim.x + threadIdx.x;
// Keep mind: Partition the input into blockDim.x subsections: i.e. for 8 threads --> 8 subsections
// collaborative load in a coalesced manner
for (int j = 0; j < SECTION_SIZE; j += blockDim.x) {
XY[threadIdx.x + j] = X[threadIdx.x + j];
}
__syncthreads();
// PHASE 1: scan inner own subsection
// At the end of this phase the last element of each subsection contains the sum of all alements in own subsection
for (int j = 1; j < SUBSECTION_SIZE; j++) {
XY[threadIdx.x * (SUBSECTION_SIZE)+j] += XY[threadIdx.x * (SUBSECTION_SIZE)+j - 1];
}
__syncthreads();
// PHASE 2: perform iterative kogge_stone_scan of the last elements of each subsections of XY loaded first in AUS
AUS[threadIdx.x] = XY[threadIdx.x * (SUBSECTION_SIZE)+(SUBSECTION_SIZE)-1];
float in;
for (unsigned int stride = 1; stride < blockDim.x; stride *= 2) {
__syncthreads();
if (threadIdx.x >= stride) {
in = AUS[threadIdx.x - stride];
}
__syncthreads();
if (threadIdx.x >= stride) {
AUS[threadIdx.x] += in;
}
}
__syncthreads();
// PHASE 3: each thread adds to its elements the new value of the last element of its predecessor's section
if (threadIdx.x > 0) {
for (unsigned int stride = 0; stride < (SUBSECTION_SIZE); stride++) {
XY[threadIdx.x * (SUBSECTION_SIZE)+stride] += AUS[threadIdx.x - 1]; // <--
}
}
__syncthreads();
// store the result into output vector
for (int j = 0; j < SECTION_SIZE; j += blockDim.x) {
Y[threadIdx.x + j] = XY[threadIdx.x + j];
}
}
如果我在块中使用一个线程和 8192 个元素,它会完美地工作,但如果我使用多个线程,则 XY[5793](或 X[5793] 完成并存储时的结果是错误的结果)。 它可以使用 4096 个元素和一个或多个线程(最多 1024 个线程)。 如果我使用 int 而不是 float 数字,它甚至可以使用一个或多个线程的 8192 个元素。
我也尝试在 MATLAB 中进行验证,这些是输出比较:
- X[5973] = 16788115 ---- MATLAB
- X[5973] = 16788114 ---- CPU
- X[5973] = 16788116 ---- GPU
正如我们所见,CPU 结果也与 MATLAB 不同,所以根据这些结果我认为问题出在浮点加法上,但我告诉你我用有序 "x.00" 浮点数(例如 {1.00, 2.00, 3.00, 4.00 ..... 8192.00})。
另一个问题是关于性能,主机代码总是比内核代码快,有这些配置参数和这些输入,正常吗?
如果你需要完整的源代码,你可以找到它here
第一次扫描可能有问题:
XY[threadIdx.x * (SUBSECTION_SIZE)+j] += XY[threadIdx.x * (SUBSECTION_SIZE)+j - 1];
这可能导致对共享内存中元素的读取不连贯。当您读取前一个元素时,不能保证任何其他线程都没有更新该值。
尝试通过将值存储在寄存器中来将此部分分成几部分。示例:
int t = XY[threadIdx.x * (SUBSECTION_SIZE)+j - 1];
__syncthreads();
XY[threadIdx.x * (SUBSECTION_SIZE)+j] += t;
8192 是 2^13
sum(1..8192) 接近8192^2/2: 8192*8193/2,比2^25多一点。
因此您需要 26 位来表示它(请参阅下面的注释)。
单精度 IEEE 754 浮点数只有 24 位有效数,因此,取决于求和的执行方式(以何种顺序),以及最终的舍入方向(通常它是默认舍入到最近的,绑定到偶数),然后结果可能不同。
注意严格来说,精确的和可以用float来表示,不用四舍五入,因为最后12位是0,所以尾数只跨越14位。但部分和不是这样。