通过基本函数在较少的大数组中聚合许多小数组
Posted
技术标签:
【中文标题】通过基本函数在较少的大数组中聚合许多小数组【英文标题】：Aggregate many small arrays in fewer large arrays by basic function 【发布时间】：2021-06-09 22:08:34 【问题描述】：我有许多小型 2D 阵列（例如 M x 32 x 40）和更少的大型 2D 阵列（例如 N x 200 x 300）。 我想将较小的矩阵“放置”在较大数组中的索引 n、i、j 处（批处理索引 n 处数组的左上索引）。这些小数组可能重叠，应该由关联和交换的函数聚合，比如加、乘等。
我认为这是很多人应该遇到的一个非常基本的场景，对吧？是否有一个 cuda 实现可以有效地支持这一点？
典型值 M = 10^6, N = 10^4
【问题讨论】：
减号运算符既不是关联也不是可交换的，因此在 CUDA 中有效地实现这一点可能比仅使用像加号运算符这样的关联+交换运算符要复杂得多。同样的事情也适用于组合运算符。这也意味着矩阵的顺序很重要。 非常好的评论！编辑了我的问题。它实际上是关于加号 不清楚你如何确定n；我认为找到一个高性能的解决方案很重要。您的每个M 数组是否都有一个与之关联的n，也许在其他数组中，它可以识别应该添加到哪个N 数组？
small n 是已知的并且是已知的（i 和 j 也是如此）。是的，每个M 数组都有一个与之关联的n。
这个问题看起来类似于有限元法装配过程。如果您搜索，您应该会找到从 CUDA 早期到现在的一系列论文。我怀疑你会找到现成的库实现
【参考方案1】：
这是一个归约操作。
除了在 cmets 中表示的内容外，我将假设 M 矩阵的分布相对于它们属于哪个 N 矩阵而言是相对均匀的，即均匀分布。这意味着对于给定的维度，将有大约 100 个 M 矩阵打算更新 N 矩阵 0，100 个用于更新 N 矩阵 1，依此类推。此外，如果我们检查 n 数组，我们会观察到索引的均匀随机模式（即没有聚集或分组）。
鉴于这对我来说可能是第一次，我将建议使用来自here 的管道的锁定/临界区算法。每个线程块将采用M 数组之一，并尝试获取锁，以便它可以更新适当的N 数组。完成后，释放锁。
我也考虑了其他方法，其中一些在代码中很明显。无论如何，对于所述条件，基于锁的方法在我的 V100 GPU 上的内核运行时间约为 40 毫秒，这是我观察到的最好的。
我还要注意，所述尺寸导致数据工作集约为 8GB。这不是问题，请注意是否在您的笔记本电脑 GPU 上按原样运行此代码。
这是一个例子：
$ cat t34.cu
#include <iostream>
#include <cstdlib>
const int N = 10000;
const int M = 1000000;
const int Mx = 32;
const int My = 40;
const int Nx = 200;
const int Ny = 300;
const int nTPB = 256;
template <typename T>
__host__ __device__
T reduction_op(T &a, const T &b) return a+b;
template <typename T>
__global__ void k(const T * __restrict__ M, T * __restrict__ N, const int * __restrict__ n, const int * __restrict__ i, const int * __restrict__ j, const int num_M)
for (int ii = 0; ii < num_M; ii++)
if (n[ii] == blockIdx.x)
for (int jj = threadIdx.x; jj < Mx*My; jj += blockDim.x)
int y = jj/Mx;
int x = jj - (y*Mx);
N[blockIdx.x*Nx*Ny + i[ii] + (j[ii]+y)*Nx + x] = reduction_op(
N[blockIdx.x*Nx*Ny + i[ii] + (j[ii]+y)*Nx + x], M[ii*Mx*My + y*Mx + x]);
__syncthreads();
// assumes Ny is whole-number divisible by sl
template <typename T>
__global__ void ki(const T * __restrict__ M, T * __restrict__ N, const int * __restrict__ n, const int * __restrict__ i, const int * __restrict__ j, const int num_M, const int sl)
extern __shared__ T s[];
for (int c = 0; c < Ny; c+=sl) // process per chunk of N array
// load shared
for (int t = threadIdx.x; t < sl*Nx; t += blockDim.x) s[t] = N[blockIdx.x*Nx*Ny + c*Nx + t];
__syncthreads();
// process chunk stack
for (int ii = 0; ii < num_M; ii++) // iterate through "stack"
if ((n[ii] == blockIdx.x) && (j[ii] < (c+sl)) && ((j[ii]+My) > c))
for (int jj = threadIdx.x; jj < sl*Mx; jj += blockDim.x)
int y = jj/Mx;
int x = jj - (y*Mx);
//y += c;
if ((y+c >= j[ii]) && (y+c < (j[ii]+My)))
s[y*Nx+x+i[ii]] = reduction_op(s[y*Nx+x+i[ii]], M[ii*Mx*My + (y+c-j[ii])*Mx + x]);
__syncthreads();
// save shared
for (int t = threadIdx.x; t < sl*Nx; t += blockDim.x) N[blockIdx.x*Nx*Ny + c*Nx + t] = s[t];
template <typename T>
__global__ void ka(const T * __restrict__ M, T * __restrict__ N, const int * __restrict__ n, const int * __restrict__ i, const int * __restrict__ j, const int num_M)
int x = threadIdx.x;
for (int y = threadIdx.y; y < My; y += blockDim.y)
atomicAdd(N+n[blockIdx.x]*Nx*Ny+(j[blockIdx.x]+y)*Nx+i[blockIdx.x]+x, M[blockIdx.x*Mx*My+y*Mx+x]);
__device__ void acquire_semaphore(volatile int *lock)
while (atomicCAS((int *)lock, 0, 1) != 0);
__device__ void release_semaphore(volatile int *lock)
*lock = 0;
__threadfence();
template <typename T>
__global__ void kl(const T * __restrict__ M, T * __restrict__ N, const int * __restrict__ n, const int * __restrict__ i, const int * __restrict__ j, const int num_M, int * __restrict__ locks)
if ((threadIdx.x == 0) && (threadIdx.y == 0))
acquire_semaphore(locks+n[blockIdx.x]);
__syncthreads();
//begin critical section
int x = threadIdx.x;
for (int y = threadIdx.y; y < My; y += blockDim.y)
N[n[blockIdx.x]*Nx*Ny + i[blockIdx.x] + (j[blockIdx.x]+y)*Nx + x] = reduction_op(
N[n[blockIdx.x]*Nx*Ny + i[blockIdx.x] + (j[blockIdx.x]+y)*Nx + x], M[blockIdx.x*Mx*My + y*Mx + x]);
// end critical section
__threadfence(); // not strictly necessary for the lock, but to make any global updates in the critical section visible to other threads in the grid
__syncthreads();
if ((threadIdx.x == 0) && (threadIdx.y == 0))
release_semaphore(locks+n[blockIdx.x]);
typedef float mt;
int main()
mt *d_M, *h_M, *d_N, *h_N, *r1, *r2;
int *d_n, *h_n, *d_i, *h_i, *d_j, *h_j;
h_M = new mt[M*Mx*My];
h_N = new mt[N*Nx*Ny];
r1 = new mt[N*Nx*Ny];
r2 = new mt[N*Nx*Ny];
h_n = new int[M];
h_i = new int[M];
h_j = new int[M];
cudaMalloc(&d_M, M*Mx*My*sizeof(mt));
cudaMalloc(&d_N, N*Nx*Ny*sizeof(mt));
cudaMalloc(&d_n, M*sizeof(int));
cudaMalloc(&d_i, M*sizeof(int));
cudaMalloc(&d_j, M*sizeof(int));
for (int i = 0; i < M; i++)
h_n[i] = rand()%N;
h_i[i] = rand()%(Nx - Mx);
h_j[i] = rand()%(Ny - My);
for (int i = 0; i < N*Nx*Ny; i++) h_N[i] = (mt)(i%3);
for (int i = 0; i < M*Mx*My; i++) h_M[i] = (mt)((i%3)+1);
cudaMemcpy(d_M, h_M, M*Mx*My*sizeof(mt), cudaMemcpyHostToDevice);
cudaMemcpy(d_N, h_N, N*Nx*Ny*sizeof(mt), cudaMemcpyHostToDevice);
cudaMemcpy(d_n, h_n, M*sizeof(int), cudaMemcpyHostToDevice);
cudaMemcpy(d_i, h_i, M*sizeof(int), cudaMemcpyHostToDevice);
cudaMemcpy(d_j, h_j, M*sizeof(int), cudaMemcpyHostToDevice);
#ifdef USE_SINGLE_N
cudaMemset(d_n, 0, M*sizeof(int));
#endif
#if 0
const int sl = 40;
const int sb = sl * Nx * sizeof(mt);
ki<<<N, nTPB, sb>>>(d_M, d_N, d_n, d_i, d_j, M, sl);
cudaMemcpy(r2, d_N, N*Nx*Ny*sizeof(mt), cudaMemcpyDeviceToHost);
#endif
dim3 block(Mx, 8);
#if 0
ka<<<M, block>>>(d_M, d_N, d_n, d_i, d_j, M);
cudaMemcpy(r2, d_N, N*Nx*Ny*sizeof(mt), cudaMemcpyDeviceToHost);
#endif
int *d_locks;
cudaMalloc(&d_locks, N*sizeof(int));
cudaMemset(d_locks, 0, N*sizeof(int));
kl<<<M, block>>>(d_M, d_N, d_n, d_i, d_j, M, d_locks);
cudaMemcpy(r2, d_N, N*Nx*Ny*sizeof(mt), cudaMemcpyDeviceToHost);
cudaMemcpy(d_N, h_N, N*Nx*Ny*sizeof(mt), cudaMemcpyHostToDevice);
k<<<N, nTPB>>>(d_M, d_N, d_n, d_i, d_j, M);
cudaMemcpy(r1, d_N, N*Nx*Ny*sizeof(mt), cudaMemcpyDeviceToHost);
for (int i = 0; i < N*Nx*Ny; i++) if (r1[i] != r2[i]) std::cout << "mismatch at: " << i << " was: " << r2[i] << " should be: " << r1[i] << std::endl; return 0;
$ nvcc -o t34 t34.cu -O3 -lineinfo
$ nvprof ./t34
==17970== NVPROF is profiling process 17970, command: ./t34
==17970== Profiling application: ./t34
==17970== Profiling result:
Type Time(%) Time Calls Avg Min Max Name
GPU activities: 34.57% 3.09036s 2 1.54518s 1.54294s 1.54742s [CUDA memcpy DtoH]
33.18% 2.96615s 1 2.96615s 2.96615s 2.96615s void k<float>(float const *, float*, int const *, int const *, int const *, int)
31.81% 2.84401s 6 474.00ms 1.4255ms 1.27035s [CUDA memcpy HtoD]
0.45% 39.949ms 1 39.949ms 39.949ms 39.949ms void kl<float>(float const *, float*, int const *, int const *, int const *, int, int*)
0.00% 2.1120us 1 2.1120us 2.1120us 2.1120us [CUDA memset]
API calls: 96.13% 8.94558s 8 1.11820s 1.9203ms 4.51030s cudaMemcpy
3.60% 334.59ms 6 55.765ms 277.58us 330.37ms cudaMalloc
0.15% 13.752ms 8 1.7190ms 1.3268ms 2.2025ms cuDeviceTotalMem
0.11% 10.472ms 808 12.959us 172ns 728.50us cuDeviceGetAttribute
0.01% 997.81us 8 124.73us 100.93us 176.73us cuDeviceGetName
0.00% 69.047us 2 34.523us 32.349us 36.698us cudaLaunchKernel
0.00% 68.013us 1 68.013us 68.013us 68.013us cudaMemset
0.00% 46.172us 8 5.7710us 1.8940us 23.025us cuDeviceGetPCIBusId
0.00% 8.5060us 16 531ns 260ns 1.5030us cuDeviceGet
0.00% 3.7870us 8 473ns 229ns 881ns cuDeviceGetUuid
0.00% 3.3980us 3 1.1320us 610ns 2.0780us cuDeviceGetCount
$
扩展讨论：
关于性能：
这是一个内存绑定算法。因此，我们可以通过确定执行操作所需的最小内存读取和写入次数，然后除以可用内存带宽来确定内核持续时间的最佳或下限，从而估计最佳内核性能。不幸的是，最小读写次数的确定取决于M 矩阵的位置，因此如果不检查n、i 和j 矩阵，一般不容易确定。
但是我们可以寻找另一种估算方法。另一种估计方法是观察每个M 矩阵更新将需要读取 2 个值并写入一个值。如果我们用它作为我们的估计，我们会得出M*Mx*My*3*sizeof(element_of_M)/GPU_memory_bandwidth。在我的 V100（~700GB/s BW）上，内核持续时间的下限约为 20 毫秒。
关于考虑的方法：
“天真”方法，内核k：每个线程块将负责N 矩阵之一，并将遍历M 矩阵，检查n 以确定M 矩阵是否会更新分配的N 矩阵。这给出了约 3 秒的非最佳运行时间，但根据 n 的分布，似乎在性能方面基本保持不变，并且可以使用“任意”缩减操作。
尝试“最佳”方法，内核ki：每个线程块将负责N 矩阵之一，但一次只会加载该矩阵的一部分。然后它将继续通过M 矩阵更新该块，类似于k 内核。这需要更多的矩阵循环，但应该“几乎”只加载或保存每个全局内存项所需的最少次数。不过，运行时间真的很长，~40s
原子方法，内核ka：每个线程块将负责M 矩阵之一，并将自动更新相关的N 矩阵。简单。运行时间“快”约为 40 毫秒。 （原子方法可能比非均匀n 分布更快。我目睹了内核运行时间低至 8 毫秒！）然而，这并不容易推广到没有原子等价物的操作，例如乘法。
基于锁的方法，内核kl：与原子方法一样，每个线程块将负责M 矩阵之一，并将首先获取相关N 矩阵上的锁。锁意味着不需要原子。对于所呈现的均匀分布的n 情况，它的性能与原子情况大致相同。它的好处是它可以轻松处理其他归约操作，例如乘法。一个缺点是，在 n 中存在非均匀随机分布的情况下，性能可能会受到影响，最坏的情况是在幼稚内核（3-5 秒）的范围内。
总体而言，如果可以放弃对任意归约运算符的要求（例如，仅使用加法），那么原子方法可能是最好的。
【讨论】：
以上是关于通过基本函数在较少的大数组中聚合许多小数组的主要内容，如果未能解决你的问题，请参考以下文章
MongoDB Schema Design - 许多小文档还是更少的大文档？