GPU中有两种内存最常用--全局内存和共享内存。全局内存容量大、速度慢，延迟高，而共享内存是较小的片上内存(SM内部)，具有较低的延迟，和较快的访问速度。共享内存的用途一般为：1、线程块内通信的通道，2、用作全局内存的可编程缓存。
在Fermi和Kepler架构上，共享内存和一级缓存是同一个物理硬件，共有64KB的空间，默认分配48KB给共享内存，16KB给一级缓存。从Maxwell架构开始，共享内存独立出来作为一个单独的硬件，每个SM有64KB的空间，但每个线程块最多仍然只能分配和使用48KB的空间。
当每个线程块开始执行时，会分配给它一定数量的共享内存。这个共享内存的地址空间被线程块中所有的线程共享。它和所在的线程块具有相同生命周期。共享内存被SM中的所有常驻线程块划分，因此，共享内存是限制设备并行性的关键资源。一个核函数使用的共享内存越多，处于并发活跃状态的线程块就越少。
在CUDA中，共享内存可以静态或动态地分配。共享内存可以在核函数内部声明为局部的，也可以在核函数外部声明成全局的。CUDA支持一维、二维和三维共享内存数组的声明。
共享内存变量用下列修饰符进行声明：
静态声明一块共享内存的语法为：
1
|
__shared__ float cache[32][32];
|
静态声明的共享内存其长度在编译时必须已知，也就是和c++的静态数组规定一致。动态声明时不指定内存的长度，在运行时指定长度，和c++中的动态数组类似，其语法为：
1
2
3
4
5
6
|
__global__ void kernel( /* ... */ ) {
extern __shared__ float cache[];
// ...
}
kernel<<<grid, block, m * sizeof(float)>>>( /* ... */ );
|
在执行核函数时必须指定核函数内部动态分配的共享内存的字节数，它作为核函数启动语句的<<<>>>的第三个参数传入。此处所说的字节数是一个线程块所需要的字节数，不是整个线程网格需要的字节数，所以上面代码中的m一般是和block相关的一个数。注意，动态声明的共享内存只能是一维的。
共享内存声明于核函数内部时，其生命周期与核函数启动的线程块的生命周期相同，如果声明于核函数外部则与程序的生命周期相同，全局的共享内存只能静态声明。全局的共享内存很少使用，此节不再过多介绍。
局部声明的共享内存无论是动态的还是静态的，它们的生命周期都和线程块相同，线程块被调度到一相SM上时为它分配共享内存，执行结束时回收为它分配的共享内存。如果SM上剩余的共享内存数量不能满足线程块的要求，则线程块不会被调度到这个SM上。
为了获得高内存带宽，共享内存被分为32个同样大小存储体，采用低位交叉编址，因此连续32个地址位于不同的存储体内，它们的数据可以在一个内存事务的时间内取到。如果多个地址请求落在相同的内存存储体中时，就会发生存储体冲突，这会导致访问请求被重复执行。
当线程束发出共享内存访问请求时，可能有三种情况：1、并行访问：多个地址访问多个存储体；2、串行访问：多个地址访问同一个存储体；3、广播访问：单一地址读取单一存储体。并行访问是比较常见的情况，最好的情况是线程束内的32个线程访问的地址都在不同的存储体内，只需要一个内存事务的时间即可取到所有数据。最坏的情况是线程束32个地址都不相同，且位于同一个存储体，需要32个内存事务才能取到所有数据。
前两种情况都是最优情况，而最后一种如果同一存储体访问的是相同的地址则是最优情况，否则不是。
低位交叉编址是将低地址用于选存储体，将高地址用于在存储体内选数据。64KB的内存需要16位地址，每个存储体有2KB内存，需要11位地址。如果从低位到高位为第0位到第15位，那么可以这样设计地址空间：
-
完全的低位交叉编址。如果将低5(0-4)位用于选存储体，高11(5-15)位用于存储体内选数据，那么连续的1字节在不同的存储体内，显然不合适，因为取一个int或者float就需要访问4个存储体。
-
如果将第2-6位用于选存储体，其余位用于存储体内选数据，则每过4字节需要跳到下一个存储体，一个int或者float在一个存储体内，相邻的int在相邻的存储体内，很合理。
-
如果将第3-7位用于选存储体，其余位用于存储体内选数据，则每过8字节需要跳到下一个存储体，一个double也将在一个存储体内，相邻的double在相邻的存储体内，也很合理。
CUDA提供第二、三两种方案，可以通过以下函数来查询或设置：
1
2
|
cudaError_t cudaDeviceGetSharedMemConfig(cudaSharedMemConfig **pConfig);
cudaError_t cudaDeviceSetSharedMemConfig(cudaSharedMemConfig config);
|
cudaSharedMemConfig的值可以为：
cudaSharedMemBankSizeDefault: device default (currently, four bytes)
cudaSharedMemBankSizeFourByte: four bytes natively
cudaSharedMemBankSizeEightByte: eight bytes natively
从Kepler架构开始，共享内存可以在一个内存事务提供64位的带宽，因此，如果将共享内存的宽度设置为8字节，则相邻的两个int或者float虽然地址不同，又在同一个存储体内，但是可以在一个内存事务中取到，因此不会发生存储体冲突。不过设置成8字节也不总是好的，它也有可能使原本位于不同存储体的两个int变成在同一个存储体中，从而发生冲突，比如
1
|
__shared__ int cache[1024];
|
在宽度为4字节时(下标对32取模)，cache[0]在0号存储上，cache[65]在1号存储体上，而在宽度为8字节时(下标除以2再对32取模)，它们都在0号存储体上，它们的访问将发生冲突。因此是否调整共享内存的宽度应当看核函数执行时各线程的访存位置。
除了配置合适的共享内存宽度以外，避免存储体冲突还可以在数组内填入一些空白位来避免冲突。
共享内存在线程块内共享，是线程块内的重要通信手段，不同线程对同一块内存的读写将引起争用。对同一块内存的读写，包括共享内存、全局内存等，CUDA采取弱一致性模型，也就是只保证同时读写的线程不会读到意料之外的值，比如一块内存本来为3，一个线程写-7，另一个线程要么读到3，要么读到-7，而不会读到它们二进制位混合的结果，比如-3，但弱一致性模型不保证读写的顺序，也就是一个线程的写入顺序对其他线程可见时，它可能和写操作被执行的实际顺序不一致。因此多个线程对同一个内存地址的读写要考虑使用原子操作进行互斥，或者使用同步操作进行同步。此节仅介绍同步操作。
同步操作用于强制程序以一个确定的顺序执行，一般为障碍，核函数中可以使用的障碍有：
1
2
3
4
|
void __syncthreads();
void __threadfence_block();
void __threadfence();
void __threadfence_system();
|
第一个操作为线程块障碍，用于将同一个线程块内的线程阻塞，直到此块内所有线程都调用此函数。注意，如果__syncthreads()出现在条件分支中，如果线程块内不是所有线程都会通过此分支，就会导致核函数无限期地阻塞，无法继续向下推进，比如归约中的同步操作误写在if分支内部：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
|
__global__ void ReductionKernel(float *vec, float *partial, int N) {
int tid = threadIdx.x + blockIdx.x * blockDim.x;
// 首先每个元素计算部分元素加和
float t = 0;
for (int i = tid; i < N; i += gridDim.x * blockDim.x)
t += vec[i];
vec[tid] = t;
__syncthreads();
// 然后在线程块内进行并行归约
if (blockDim.x >= 1024 && tid < 512) {
vec[tid] += vec[tid + 512];
__syncthreads(); // 核函数阻塞于此
}
if (blockDim.x >= 512 && tid < 256) {
vec[tid] += vec[tid + 256];
__syncthreads();
}
// ...
|
__threadfence_block等函数是内存栅栏，用于确保栅栏前的任何内存读写操作对栅栏后的线程都是可见的。__threadfence_block是线程块级的栅栏，保存块内线程对内存写操作的同步，__threadfence是线程网格级的栅栏，__threadfence是系统级的栅栏，系统级的栅栏可以保证多个设备对主机内存的读写的同步。一般而言，最好不好使用内存栅栏，它可能引起很大的性能损失。
另一个需要注意的地方是线程对内存的更改可能不知情，比如：
1
2
3
4
5
6
|
int global_num = 1;
int f1() {
global_num = 2;
std::cout << global_num << std::endl;
}
|
f1对global_num赋值后马上读取，在单线程下，这个值显然是刚刚写入的2，因此没有必要再去读一次内存，编译器也是这样认为的，极有可能会将输出语句中对global_num的读取操作优化掉。但是在多线程的情况下，global_num可能被其它线程更改，因此编译器的优化可能是错误的。为此，c、c++和CUDA都引入了关键字volatile来表示一块内存可能因为其它原因而更改，指示编译器在优化时不要对它作任何假设，这个关键字在《CUDA教程3》一节的归约程序，和《CUDA教程4》一节的矩阵向量乘法中已经使用过了。
《CUDA教程3》一节给出了一个归约程序：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
|
// 需要假定N足够大，否则靠后的线程块内对vec的访问可能越界
__global__ void ReductionKernel(float *vec, float *partial, int N) {
int tid = threadIdx.x + blockIdx.x * blockDim.x;
// 首先每个元素计算部分元素加和
float t = 0;
for (int i = tid; i < N; i += gridDim.x * blockDim.x)
t += vec[i];
vec[tid] = t;
__syncthreads();
// 然后在线程块内进行并行归约
if (blockDim.x >= 1024 && tid < 512)
vec[tid] += vec[tid + 512];
__syncthreads();
if (blockDim.x >= 512 && tid < 256)
vec[tid] += vec[tid + 256];
__syncthreads();
if (blockDim.x >= 256 && tid < 128)
vec[tid] += vec[tid + 128];
__syncthreads();
if (blockDim.x >= 128 && tid < 64)
vec[tid] += vec[tid + 64];
__syncthreads();
if (blockDim.x >= 64 && tid < 32)
vec[tid] += vec[tid + 32];
__syncthreads();
if (tid < 32) {
volatile float *vmem = vec;
vmem[tid] += vmem[tid + 32];
vmem[tid] += vmem[tid + 16];
vmem[tid] += vmem[tid + 8];
vmem[tid] += vmem[tid + 4];
vmem[tid] += vmem[tid + 2];
vmem[tid] += vmem[tid + 1];
}
// 0号线程将块内求和的结果写出
if (threadIdx.x == 0)
partial[blockIdx.x] = vec[tid];
}
|
这个程序有两个问题：1、如果输入数组的长度比较短，那么vec[tid]可能越界；2、输入数组会被更改。第1个问题可以通过添加一些条件判断来解决，比如是每次写vec时判断下标是否小于N。而第二个问题就比较严重了，比如计算一个向量的2范数，许多时候我们都不想破坏这个向量，那么只好将它备份一下了，这势必造成内存和时间的浪费，而且对全局内存的访问比较慢，延迟比较高。显然使用共享内存是一个绝佳的解决方案。大致代码为：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
|
__global__ void ReductionKernel(const float *vec, float *partial, int N) {
extern __shared__ float cache[];
int tid = threadIdx.x + blockIdx.x * blockDim.x;
// 首先每个元素计算部分元素加和
float t = 0;
for (int i = tid; i < N; i += gridDim.x * blockDim.x)
t += vec[i];
cache[threadIdx.x] = t; // 不能再用tid了，共享内存只是线程块内共享，线程块间不是同一个
__syncthreads();
// 然后在线程块内进行并行归约
if (blockDim.x >= 1024 && tid < 512)
cache[threadIdx.x] += cache[threadIdx.x + 512];
__syncthreads();
if (blockDim.x >= 512 && tid < 256)
cache[threadIdx.x] += cache[threadIdx.x + 256];
__syncthreads();
if (blockDim.x >= 256 && tid < 128)
cache[threadIdx.x] += cache[threadIdx.x + 128];
__syncthreads();
if (blockDim.x >= 128 && tid < 64)
cache[threadIdx.x] += cache[threadIdx.x + 64];
__syncthreads();
if (blockDim.x >= 64 && tid < 32)
cache[threadIdx.x] += cache[threadIdx.x + 32];
__syncthreads();
// 线程束内归约
if (tid < 32) {
volatile float *vmem = cache;
vmem[threadIdx.x] += vmem[threadIdx.x + 32];
vmem[threadIdx.x] += vmem[threadIdx.x + 16];
vmem[threadIdx.x] += vmem[threadIdx.x + 8];
vmem[threadIdx.x] += vmem[threadIdx.x + 4];
vmem[threadIdx.x] += vmem[threadIdx.x + 2];
vmem[threadIdx.x] += vmem[threadIdx.x + 1];
}
// 0号线程将块内求和的结果写出
if (threadIdx.x == 0)
partial[blockIdx.x] = cache[0];
}
|
在线程束内归约之前，每个线程本来就没有读取过cache对应位置附近的值，编译器不可能将读取内存的指令优化掉，但线程束内归约就不同了，编译器可能优化掉部分内存读取指令，导致线程读不到新值而产生错误的计算结果。最后将线程块归约的结果写出时也可能会优化掉内存读取指令，但即使优化掉，结果也是正确的，所以不需要volatile修饰。
上面的核函数使用动态分配的共享内存，调用的代码大致为：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
|
int main() {
std::vector<float> nums(N); // 保存需要求和的数据
int block_size = 1024; // 线程块内有1024个线程
int grid_size = 26; // k20m有13个sm，每个sm可容纳2048个线程，总共26个线程块
float *vec = nullptr;
float *partial = nullptr; // 保存每个线程块的求和结果
cudaMalloc(&vec, sizeof(float) * N);
cudaMalloc(&partial, sizeof(float) * grid_size);
cudaMemcpy(vec, nums.data(), sizeof(float) * N, cudaMemcpyHostToDevice);
// 线程块内每个线程一个共享内存单元，总共block_size个float
ReductionKernel<<<grid_size, block_size, block_size * sizeof(float)>>>
(dvec, dpartial, N);
std::vector<float> temp(grid_size);
cudaMemcpy(temp.data(), partial, sizeof(float) * grid_size,
cudaMemcpyDeviceToHost);
float sum = 0;
for (int i = 0; i < grid_size; ++i)
sum += partial[i];
cudaFree(vec);
cudaFree(partial);
return 0;
}
|
我们也可以使用模板进行优化：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
|
template<int BlockSize>
__global__ void ReductionKernel(const float *vec, float *partial, int N) {
__shared__ float cache[BlockSize];
int tid = threadIdx.x + blockIdx.x * blockDim.x;
// 首先每个元素计算部分元素加和
float t = 0;
for (int i = tid; i < N; i += gridDim.x * blockDim.x)
t += vec[i];
cache[threadIdx.x] = t; // 不能再用tid了，共享内存只是线程块内共享，线程块间不是同一个
__syncthreads();
// 然后在线程块内进行并行归约
if (BlockSize >= 1024 && tid < 512)
cache[threadIdx.x] += cache[threadIdx.x + 512];
__syncthreads();
if (BlockSize >= 512 && tid < 256)
cache[threadIdx.x] += cache[threadIdx.x + 256];
__syncthreads();
if (BlockSize >= 256 && tid < 128)
cache[threadIdx.x] += cache[threadIdx.x + 128];
__syncthreads();
if (BlockSize >= 128 && tid < 64)
cache[threadIdx.x] += cache[threadIdx.x + 64];
__syncthreads();
if (BlockSize >= 64 && tid < 32)
cache[threadIdx.x] += cache[threadIdx.x + 32];
__syncthreads();
// 线程束内归约
if (tid < 32) {
volatile float *vmem = cache;
vmem[threadIdx.x] += vmem[threadIdx.x + 32];
vmem[threadIdx.x] += vmem[threadIdx.x + 16];
vmem[threadIdx.x] += vmem[threadIdx.x + 8];
vmem[threadIdx.x] += vmem[threadIdx.x + 4];
vmem[threadIdx.x] += vmem[threadIdx.x + 2];
vmem[threadIdx.x] += vmem[threadIdx.x + 1];
}
// 0号线程将块内求和的结果写出
if (threadIdx.x == 0)
partial[blockIdx.x] = cache[0];
}
constexpr int block_size = 512;
ReductionKernel<block_size><<<grid_size, block_size>>>(dvec, dpartial, N);
|
这样的好处在于编译器生成代码时可以根据block_size将不需要的if语句优化掉，缺点是模板参数必须是编译器常量，即编译时就必须确定线程块的大小。
《CUDA教程4》一节中有一个线程束内归约的矩阵向量乘法，代码大致为：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
|
__global__ void matrixVectorMulKernel(
double *out_vec,
volatile double *cache,
const double *mat,
const double *in_vec,
const int M,
const int N
) {
unsigned int tid = blockIdx.x * blockDim.x + threadIdx.x;
int row = tid / 32;
if (row < M) {
double t = 0;
for (int col = tid % 32; col < N; col += 32)
t += mat[row * N + col] * in_vec[col];
cache[tid] = t;
if (tid % 32 < 16) {
cache[tid] += cache[tid + 16];
cache[tid] += cache[tid + 8];
cache[tid] += cache[tid + 4];
cache[tid] += cache[tid + 2];
cache[tid] += cache[tid + 1];
}
if (tid % 32 == 0)
out_vec[row] = cache[tid];
}
}
int main() {
double *mat = nullptr;
double *vec = nullptr;
double *result = nullptr;
cudaMalloc(&mat, M * N * sizeof(double));
cudaMalloc(&vec, N * sizeof(double));
cudaMalloc(&result, M * sizeof(double));
// generate mat and vec
dim3 block{512};
dim3 grid{(M * 32 + block.x - 1) / block.x};
double *cache = nullptr;
cudaMalloc(&result, block.x * grid.x * sizeof(double));
MatMulVecKernel<<<grid, block>>>(result, mat, vec, M, N);
// use result
cudaFree(mat);
cudaFree(vec);
cudaFree(result);
}
|
这段代码中，我们额外分配了一段内存cache用来保存归约过程中的数据，显然将它替换成共享内存也是一个非常好的改进：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
|
template<int BlockSize>
__global__ void matrixVectorMulKernel(
double *out_vec,
const double *mat,
const double *in_vec,
const int M,
const int N
) {
volatile __shared__ double cache[BlockSize];
unsigned int tid = blockIdx.x * blockDim.x + threadIdx.x;
int row = tid / 32;
if (row < M) {
double t = 0;
for (int col = tid % 32; col < N; col += 32)
t += mat[row * N + col] * in_vec[col];
cache[threadIdx.x] = t; // 不能用tid了
if (tid % 32 < 16) {
cache[threadIdx.x] += cache[threadIdx.x + 16];
cache[threadIdx.x] += cache[threadIdx.x + 8];
cache[threadIdx.x] += cache[threadIdx.x + 4];
cache[threadIdx.x] += cache[threadIdx.x + 2];
cache[threadIdx.x] += cache[threadIdx.x + 1];
}
if (tid % 32 == 0)
out_vec[row] = cache[threadIdx.x];
}
}
int main() {
double *mat = nullptr;
double *vec = nullptr;
double *result = nullptr;
cudaMalloc(&mat, M * N * sizeof(double));
cudaMalloc(&vec, N * sizeof(double));
cudaMalloc(&result, M * sizeof(double));
// generate mat and vec
constexpr int block_size = 512;
dim3 block{block_size};
dim3 grid{(M * 32 + block.x - 1) / block.x};
MatMulVecKernel<block_size><<<grid, block>>>(result, mat, vec, M, N);
// use result
cudaFree(mat);
cudaFree(vec);
cudaFree(result);
}
|
对于线程束内的归约，使用线程块级的通信设施总是显得很浪费。CUDA从Kepler架构引入了一组指令用于线程束内通信，称为线程束洗牌指令(warp shuffle functions)。线程束洗牌指令有四种，在CUDA9.0以前，它们为：
1
2
3
4
5
6
7
8
9
10
11
|
int __shfl(int var, int srcLane, int width = warpSize);
float __shfl(float var, int srcLane, int width = warpSize);
int __shfl_up(int var, unsigned int delta, int width = warpSize);
float __shfl_up(float var, unsigned int delta, int width = warpSize);
int __shfl_down(int var, unsigned int delta, int width = warpSize);
float __shfl_down(float var, unsigned int delta, int width = warpSize);
int __shfl_xor(int var, int laneMask, int width = warpSize);
float __shfl_xor(float var, int laneMask, int width = warpSize);
|
文档上列出了这8个函数，但实际上double版本也是可用的。这些函数在CUDA9.0中标记为deprecated，在CUDA10中移除，取而代之的是：
1
2
3
4
|
T __shfl_sync(unsigned mask, T var, int srcLane, int width = warpSize);
T __shfl_up_sync(unsigned mask, T var, unsigned int delta, int width = warpSize);
T __shfl_down_sync(unsigned mask, T var, unsigned int delta, int width = warpSize);
T __shfl_xor_sync(unsigned mask, T var, int laneMask, int width = warpSize);
|
其中T可以是int, unsigned int, long, unsigned long, long long, unsigned long long, float和double，如果包含了头文件cuda_fp16.h，T还可以是__half和__half2。多出来的参数mask用于标记线程束内哪些线程为活跃线程，当它为0xffffffff时所有线程都是活跃线程。
__shfl用于将线程束内第srcLane个线程中的var的值广播到线程束内其它线程。__shfl_up用于将值传递给线程id为自己的id+delta的线程：
__shfl_down用于将值传递给线程id为自己的id-delta的线程：
__shlf_xor用于交换值：
当laneMask取不同的值时，效果大概为：
mask: 1
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18...
1 0 3 2 5 4 7 6 9 8 11 10 13 12 15 14 17 16 19...
mask: 2
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18...
2 3 0 1 6 7 4 5 10 11 8 9 14 15 12 13 18 19 16...
mask: 4
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18...
4 5 6 7 0 1 2 3 12 13 14 15 8 9 10 11 20 21 22...
mask: 8
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18...
8 9 10 11 12 13 14 15 0 1 2 3 4 5 6 7 24 25 26...
mask: 16
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18...
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 0 1 2 ...
例如使用线程束洗牌指令实现矩阵向量乘法：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
|
__global__ void matrixVectorMulKernel(
double *out_vec,
const double *mat,
const double *in_vec,
const int M,
const int N
) {
unsigned int tid = blockIdx.x * blockDim.x + threadIdx.x;
int row = tid / 32;
if (row < M) {
double t = 0;
for (int col = tid % 32; col < N; col += 32)
t += mat[row * N + col] * in_vec[col];
t += __shfl_xor_sync(0xffffffff, t, 16);
t += __shfl_xor_sync(0xffffffff, t, 8);
t += __shfl_xor_sync(0xffffffff, t, 4);
t += __shfl_xor_sync(0xffffffff, t, 2);
t += __shfl_xor_sync(0xffffffff, t, 1);
if (tid % 32 == 0)
out_vec[row] = t;
}
}
int main() {
double *mat = nullptr;
double *vec = nullptr;
double *result = nullptr;
cudaMalloc(&mat, M * N * sizeof(double));
cudaMalloc(&vec, N * sizeof(double));
cudaMalloc(&result, M * sizeof(double));
// generate mat and vec
dim3 block{512};
dim3 grid{(M * 32 + block.x - 1) / block.x};
MatMulVecKernel<<<grid, block>>>(result, mat, vec, M, N);
// use result
cudaFree(mat);
cudaFree(vec);
cudaFree(result);
}
|
非常漂亮，不需要共享内存，不需要模板，性能也非常好。
1、如何实现矩阵转置？
1
2
3
4
5
6
7
8
9
10
11
12
|
__global__ void transposeKernelNaive(
float *At,
const float *A,
const int M,
const int N
) {
unsigned int ix = threadIdx.x + blockIdx.x * blockDim.x;
unsigned int iy = threadIdx.y + blockIdx.y * blockDim.y;
if (ix < M && iy < N) {
At[ix * N + iy] = A[iy * M + ix];
}
}
|
对全局内存的访问，读可以合并，但写不合并，如何改进？
2、使用线程束洗牌指令做线程束内归约还可以怎么做，本节最后使用线程束洗牌指令实现矩阵向量乘法的程序是否还可以改进？
3、进阶：矩阵乘法如何实现？
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
|
__global__ void MatMulKernelNaive(
float *C, // M x N
const float *A, // M x K
const float *B, // K x N
const int M,
const int N,
const int K
) {
unsigned int ix = blockIdx.x * blockDim.x + threadIdx.x;
unsigned int iy = blockIdx.y * blockDim.y + threadIdx.y;
if (ix < M && iy < N) {
float t = 0;
for (int k = 0; k < K; ++k)
t += A[ix * K + k] * B[k * N + iy];
C[ix * N + iy] = t;
}
}
|
对A、B的访问都不能合并，如何改进？