CUDA：带有 3D 内核的嵌套 FOR 循环：如何确定线程应写入结果的位置？
CUDA: nested FOR-loops with 3D kernel: How to determine the position where threads should write the result?
在接下来的迭代中，向量“result”中的位置由“counter”决定。请同时查看 FOR 循环中的起始值，这表明遗漏了一些迭代。
int counter = (int)0;
for (z=0; z<N; z++)
for (y=z; y<N; y++)
for (x=y; x<N; x++)
{
result[counter] = A[z] + A[y] + A[x];
counter++;
}
但是，如果我将此迭代转换为 3d 内核，我想知道如何为每个线程提供正确的位置以写入向量“结果”？是否有计算出缺失线程的数学解决方案？
int idz = blockIdx.z*blockDim.z + threadIdx.z;
int idy = blockIdx.y*blockDim.y + threadIdx.y;
int idx = blockIdx.x*blockDim.x + threadIdx.x;
if(idx >= idy && idy >= idz)
{
result[?????] = A[idz] + A[idy] + A[idx];
}
我试图用高斯方程跳过缺失值：
0.5 * idx * (idx + 1)
但我离解决方案还很远。有人有想法吗？
我想解决这个问题的方法可能不止一种。我将概述一种方法。
前n个整数之和和前n个整数之和的平方和的公式为readily available为：
Σi i:1..n = (n)(n+1)/2
Σi^2 i:1..n = (n)(n+1)(2n+1)/6
让我们从构建一个 table 开始，显示 x 中每个 z,y 位置的循环迭代次数。如果我们为 N 选择一个小数字，例如 4:
，那么将其可视化会有所帮助
z| 0 1 2 3
y
-
0 4
1 3 3
2 2 2 2
3 1 1 1 1
从上面的table可以看出，任意一列的和都是前n个整数的和，其中n可以从N和z索引推导出来. (n = N-z).
我们可以使用这个观察结果，加上前 n 个整数之和的公式（如上所示），来计算给定 z 值的偏移量（到结果数组中）。此偏移量将是 z 左侧列的总和。从上面的模式我们可以观察到z索引贡献的offset是：
Σ(i)(i+1)/2 i: (N-z)..N
即：
1/2 * ( Σ(i^2) + Σi) i: (N-z)..N
现在，对 (N-z)..N 求和很不方便，因此要使用上面给出的公式来实现这一点，我们应将每个求和重写为两个求和的差值，一个从 1..N 减去一个来自 1..z:
1/2 * ( Σ(i^2)|1..N - Σ(i^2)|1..z + Σi|1..N - Σi|1..z) )
如果我们将开头给出的总和的公式代入上述表达式，我们可以得出 z 索引贡献的偏移量的封闭式表达式。
(((N*(N+1)*(2*N + 1) - (N-z)*(N-z+1)*(2*(N-z)+1))/3 + (N*(N+1) - (N-z)*(N-z+1)))/2
现在我们需要为 y 索引贡献的偏移量（进入结果数组）提出一个公式。同样，通过检查前面介绍的 z,y 数组中的模式，我们可以观察到这是由以下公式给出的：
Σi i: (N-y)..(N-z)
我们同样可以将其转换为两个总和的差值进行计算，然后使用我们的公式计算整数和：
((N-z)(N-z+1) - (N-y)(N-y+1))/2
我们还必须对 x 索引进行调整。
这是一个完整的例子：
$ cat t112.cu
#include <iostream>
template <typename T>
void on_cpu(int N, T *A, T *result){
int counter = 0;
int x,y,z;
for (z=0; z<N; z++)
for (y=z; y<N; y++)
for (x=y; x<N; x++){
result[counter] = A[z] + A[y] + A[x];
counter++;}
std::cout << " counter = " << counter << std::endl;
}
__host__ __device__ int compute_z_offset(int N, int z){
int i = N-z;
// 0.5*(sum n^2|1..N - sum n^2|1..i + sum n|1..N - sum n|1..i )
int sumn2 = ((N*(N+1)*(2*N + 1)) - (i*(i+1)*(2*i + 1)))/3;
int sumn = ((N*(N+1)) - (i*(i+1)));
return (sumn2 + sumn)>>2;
}
__host__ __device__ int compute_y_offset(int N, int z, int y){
int i = N-z;
int j = N-y;
return (i*(i+1) - j*(j+1))>>1;
}
template <typename T>
__global__ void on_gpu(int N, T *A, T *result){
int idz = blockIdx.z*blockDim.z + threadIdx.z;
int idy = blockIdx.y*blockDim.y + threadIdx.y;
int idx = blockIdx.x*blockDim.x + threadIdx.x;
if(idx >= idy && idy >= idz && idx < N && idy < N && idz < N){
result[compute_z_offset(N, idz) + compute_y_offset(N, idz, idy) + (idx-idy)] = A[idz] + A[idy] + A[idx];
}
}
int main(){
const int N = 9;
int *A, *result, *d_A, *d_result, *h_result;
A = new int[N];
result = new int[N*N*N];
h_result = new int[N*N*N];
cudaMalloc(&d_result, N*N*N*sizeof(result[0]));
cudaMalloc(&d_A, N*sizeof(A[0]));
for (int i = 0; i < N; i++) A[i] = i+1;
cudaMemcpy(d_A, A, N*sizeof(A[0]), cudaMemcpyHostToDevice);
cudaMemset(d_result, 0, N*N*N*sizeof(result[0]));
on_cpu(N, A, result);
dim3 block(8,8,8);
dim3 grid((N+block.x-1)/block.x, (N+block.y-1)/block.y, (N+block.z-1)/block.z);
on_gpu<<<grid, block>>>(N, d_A, d_result);
cudaMemcpy(h_result, d_result, N*N*N*sizeof(result[0]), cudaMemcpyDeviceToHost);
for (int i = 0; i < N*N*N; i++) if (result[i] != h_result[i]) {std::cout << "mismatch at: " << i << " was: " << h_result[i] << " should be: " << result[i] << std::endl; break;}
}
$ nvcc -o t112 t112.cu
$ cuda-memcheck ./t112
========= CUDA-MEMCHECK
counter = 165
========= ERROR SUMMARY: 0 errors
$
我还没有彻底测试过这个。我的目的是演示如何完成，但可能存在缺陷。此外，我想索引计算可能会被简化（组合类似项等）。为了方便起见，我只是将结果数组的大小设置为 N*N*N，尽管并非所有这些位置都会被“使用”。您可以使用我已经在此处提供的公式分析确定必要的大小。
这里是代码的略微编辑版本。大多数情况下，它只是被修改为允许测试各种 N 值，但我也在 z 偏移计算中做了一些类似的术语组合。它似乎可以达到 N 值 1000，这是我想要测试的范围。我还演示了如何预先计算输出数组所需的大小。
#include <iostream>
#include <cstdlib>
template <typename T>
int on_cpu(int N, T *A, T *result){
int counter = 0;
int x,y,z;
for (z=0; z<N; z++)
for (y=z; y<N; y++)
for (x=y; x<N; x++){
result[counter] = A[z] + A[y] + A[x];
counter++;}
std::cout << " N: " << N << " counter = " << counter << std::endl;
return counter;
}
__host__ __device__ int compute_z_offset(int N, int z){
int i = N-z;
// 0.5*(sum n^2|1..N - sum n^2|1..i + sum n|1..N - sum n|1..i
//int sumn2 = ((N*(N+1)*(2*N + 1)) - (i*(i+1)*(2*i + 1)))/3;
//int sumn = ((N*(N+1)) - (i*(i+1)));
//return (sumn2 + sumn)>>2;
return ((N*(N+1)*(2*N + 4)) - (i*(i+1)*(2*i + 4)))/12;
}
__host__ __device__ int compute_y_offset(int N, int z, int y){
int i = N-z;
int j = N-y;
return (i*(i+1) - j*(j+1))>>1;
}
template <typename T>
__global__ void on_gpu(int N, T *A, T *result){
int idz = blockIdx.z*blockDim.z + threadIdx.z;
int idy = blockIdx.y*blockDim.y + threadIdx.y;
int idx = blockIdx.x*blockDim.x + threadIdx.x;
if(idx >= idy && idy >= idz && idx < N && idy < N && idz < N)
{
result[compute_z_offset(N, idz) + compute_y_offset(N, idz, idy) + (idx-idy)] = A[idz] + A[idy] + A[idx];
}
}
int main(int argc, char *argv[]){
int N = 32;
if (argc > 1) N = atoi(argv[1]);
int *A, *result, *d_A, *d_result, *h_result;
A = new int[N];
int max_N_size = compute_z_offset(N, N-1) + 1;
result = new int[max_N_size];
h_result = new int[max_N_size];
cudaMalloc(&d_result, sizeof(result[0])*max_N_size);
cudaMalloc(&d_A, N*sizeof(A[0]));
for (int my_N = 1; my_N<N; my_N++){
for (int i = 0; i < my_N; i++) A[i] = i+1;
int N_size = on_cpu(my_N, A, result);
int test_N_size = compute_z_offset(my_N, my_N-1) + 1;
if (test_N_size != N_size) std::cout << "mismatch on size calc" << std::endl;
cudaMemcpy(d_A, A, sizeof(A[0])*my_N, cudaMemcpyHostToDevice);
cudaMemset(d_result, 0, sizeof(result[0])*N_size);
dim3 block(8,8,8);
dim3 grid((my_N+block.x-1)/block.x, (my_N+block.y-1)/block.y, (my_N+block.z-1)/block.z);
on_gpu<<<grid, block>>>(my_N, d_A, d_result);
cudaMemcpy(h_result, d_result, sizeof(result[0])*N_size, cudaMemcpyDeviceToHost);
for (int i = 0; i < N_size; i++) if (result[i] != h_result[i]) {std::cout << "mismatch at: " << i << " was: " << h_result[i] << " should be: " << result[i] << std::endl; break;}
}
}
在接下来的迭代中，向量“result”中的位置由“counter”决定。请同时查看 FOR 循环中的起始值，这表明遗漏了一些迭代。
int counter = (int)0;
for (z=0; z<N; z++)
for (y=z; y<N; y++)
for (x=y; x<N; x++)
{
result[counter] = A[z] + A[y] + A[x];
counter++;
}
但是，如果我将此迭代转换为 3d 内核，我想知道如何为每个线程提供正确的位置以写入向量“结果”？是否有计算出缺失线程的数学解决方案？
int idz = blockIdx.z*blockDim.z + threadIdx.z;
int idy = blockIdx.y*blockDim.y + threadIdx.y;
int idx = blockIdx.x*blockDim.x + threadIdx.x;
if(idx >= idy && idy >= idz)
{
result[?????] = A[idz] + A[idy] + A[idx];
}
我试图用高斯方程跳过缺失值：
0.5 * idx * (idx + 1)
但我离解决方案还很远。有人有想法吗？
我想解决这个问题的方法可能不止一种。我将概述一种方法。
前n个整数之和和前n个整数之和的平方和的公式为readily available为：
Σi i:1..n = (n)(n+1)/2
Σi^2 i:1..n = (n)(n+1)(2n+1)/6
让我们从构建一个 table 开始，显示 x 中每个 z,y 位置的循环迭代次数。如果我们为 N 选择一个小数字，例如 4:
z| 0 1 2 3
y
-
0 4
1 3 3
2 2 2 2
3 1 1 1 1
从上面的table可以看出，任意一列的和都是前n个整数的和，其中n可以从N和z索引推导出来. (n = N-z).
我们可以使用这个观察结果，加上前 n 个整数之和的公式（如上所示），来计算给定 z 值的偏移量（到结果数组中）。此偏移量将是 z 左侧列的总和。从上面的模式我们可以观察到z索引贡献的offset是：
Σ(i)(i+1)/2 i: (N-z)..N
即：
1/2 * ( Σ(i^2) + Σi) i: (N-z)..N
现在，对 (N-z)..N 求和很不方便，因此要使用上面给出的公式来实现这一点，我们应将每个求和重写为两个求和的差值，一个从 1..N 减去一个来自 1..z:
1/2 * ( Σ(i^2)|1..N - Σ(i^2)|1..z + Σi|1..N - Σi|1..z) )
如果我们将开头给出的总和的公式代入上述表达式，我们可以得出 z 索引贡献的偏移量的封闭式表达式。
(((N*(N+1)*(2*N + 1) - (N-z)*(N-z+1)*(2*(N-z)+1))/3 + (N*(N+1) - (N-z)*(N-z+1)))/2
现在我们需要为 y 索引贡献的偏移量（进入结果数组）提出一个公式。同样，通过检查前面介绍的 z,y 数组中的模式，我们可以观察到这是由以下公式给出的：
Σi i: (N-y)..(N-z)
我们同样可以将其转换为两个总和的差值进行计算，然后使用我们的公式计算整数和：
((N-z)(N-z+1) - (N-y)(N-y+1))/2
我们还必须对 x 索引进行调整。
这是一个完整的例子：
$ cat t112.cu
#include <iostream>
template <typename T>
void on_cpu(int N, T *A, T *result){
int counter = 0;
int x,y,z;
for (z=0; z<N; z++)
for (y=z; y<N; y++)
for (x=y; x<N; x++){
result[counter] = A[z] + A[y] + A[x];
counter++;}
std::cout << " counter = " << counter << std::endl;
}
__host__ __device__ int compute_z_offset(int N, int z){
int i = N-z;
// 0.5*(sum n^2|1..N - sum n^2|1..i + sum n|1..N - sum n|1..i )
int sumn2 = ((N*(N+1)*(2*N + 1)) - (i*(i+1)*(2*i + 1)))/3;
int sumn = ((N*(N+1)) - (i*(i+1)));
return (sumn2 + sumn)>>2;
}
__host__ __device__ int compute_y_offset(int N, int z, int y){
int i = N-z;
int j = N-y;
return (i*(i+1) - j*(j+1))>>1;
}
template <typename T>
__global__ void on_gpu(int N, T *A, T *result){
int idz = blockIdx.z*blockDim.z + threadIdx.z;
int idy = blockIdx.y*blockDim.y + threadIdx.y;
int idx = blockIdx.x*blockDim.x + threadIdx.x;
if(idx >= idy && idy >= idz && idx < N && idy < N && idz < N){
result[compute_z_offset(N, idz) + compute_y_offset(N, idz, idy) + (idx-idy)] = A[idz] + A[idy] + A[idx];
}
}
int main(){
const int N = 9;
int *A, *result, *d_A, *d_result, *h_result;
A = new int[N];
result = new int[N*N*N];
h_result = new int[N*N*N];
cudaMalloc(&d_result, N*N*N*sizeof(result[0]));
cudaMalloc(&d_A, N*sizeof(A[0]));
for (int i = 0; i < N; i++) A[i] = i+1;
cudaMemcpy(d_A, A, N*sizeof(A[0]), cudaMemcpyHostToDevice);
cudaMemset(d_result, 0, N*N*N*sizeof(result[0]));
on_cpu(N, A, result);
dim3 block(8,8,8);
dim3 grid((N+block.x-1)/block.x, (N+block.y-1)/block.y, (N+block.z-1)/block.z);
on_gpu<<<grid, block>>>(N, d_A, d_result);
cudaMemcpy(h_result, d_result, N*N*N*sizeof(result[0]), cudaMemcpyDeviceToHost);
for (int i = 0; i < N*N*N; i++) if (result[i] != h_result[i]) {std::cout << "mismatch at: " << i << " was: " << h_result[i] << " should be: " << result[i] << std::endl; break;}
}
$ nvcc -o t112 t112.cu
$ cuda-memcheck ./t112
========= CUDA-MEMCHECK
counter = 165
========= ERROR SUMMARY: 0 errors
$
我还没有彻底测试过这个。我的目的是演示如何完成，但可能存在缺陷。此外，我想索引计算可能会被简化（组合类似项等）。为了方便起见，我只是将结果数组的大小设置为 N*N*N，尽管并非所有这些位置都会被“使用”。您可以使用我已经在此处提供的公式分析确定必要的大小。
这里是代码的略微编辑版本。大多数情况下，它只是被修改为允许测试各种 N 值，但我也在 z 偏移计算中做了一些类似的术语组合。它似乎可以达到 N 值 1000，这是我想要测试的范围。我还演示了如何预先计算输出数组所需的大小。
#include <iostream>
#include <cstdlib>
template <typename T>
int on_cpu(int N, T *A, T *result){
int counter = 0;
int x,y,z;
for (z=0; z<N; z++)
for (y=z; y<N; y++)
for (x=y; x<N; x++){
result[counter] = A[z] + A[y] + A[x];
counter++;}
std::cout << " N: " << N << " counter = " << counter << std::endl;
return counter;
}
__host__ __device__ int compute_z_offset(int N, int z){
int i = N-z;
// 0.5*(sum n^2|1..N - sum n^2|1..i + sum n|1..N - sum n|1..i
//int sumn2 = ((N*(N+1)*(2*N + 1)) - (i*(i+1)*(2*i + 1)))/3;
//int sumn = ((N*(N+1)) - (i*(i+1)));
//return (sumn2 + sumn)>>2;
return ((N*(N+1)*(2*N + 4)) - (i*(i+1)*(2*i + 4)))/12;
}
__host__ __device__ int compute_y_offset(int N, int z, int y){
int i = N-z;
int j = N-y;
return (i*(i+1) - j*(j+1))>>1;
}
template <typename T>
__global__ void on_gpu(int N, T *A, T *result){
int idz = blockIdx.z*blockDim.z + threadIdx.z;
int idy = blockIdx.y*blockDim.y + threadIdx.y;
int idx = blockIdx.x*blockDim.x + threadIdx.x;
if(idx >= idy && idy >= idz && idx < N && idy < N && idz < N)
{
result[compute_z_offset(N, idz) + compute_y_offset(N, idz, idy) + (idx-idy)] = A[idz] + A[idy] + A[idx];
}
}
int main(int argc, char *argv[]){
int N = 32;
if (argc > 1) N = atoi(argv[1]);
int *A, *result, *d_A, *d_result, *h_result;
A = new int[N];
int max_N_size = compute_z_offset(N, N-1) + 1;
result = new int[max_N_size];
h_result = new int[max_N_size];
cudaMalloc(&d_result, sizeof(result[0])*max_N_size);
cudaMalloc(&d_A, N*sizeof(A[0]));
for (int my_N = 1; my_N<N; my_N++){
for (int i = 0; i < my_N; i++) A[i] = i+1;
int N_size = on_cpu(my_N, A, result);
int test_N_size = compute_z_offset(my_N, my_N-1) + 1;
if (test_N_size != N_size) std::cout << "mismatch on size calc" << std::endl;
cudaMemcpy(d_A, A, sizeof(A[0])*my_N, cudaMemcpyHostToDevice);
cudaMemset(d_result, 0, sizeof(result[0])*N_size);
dim3 block(8,8,8);
dim3 grid((my_N+block.x-1)/block.x, (my_N+block.y-1)/block.y, (my_N+block.z-1)/block.z);
on_gpu<<<grid, block>>>(my_N, d_A, d_result);
cudaMemcpy(h_result, d_result, sizeof(result[0])*N_size, cudaMemcpyDeviceToHost);
for (int i = 0; i < N_size; i++) if (result[i] != h_result[i]) {std::cout << "mismatch at: " << i << " was: " << h_result[i] << " should be: " << result[i] << std::endl; break;}
}
}