内存的管理和访问是现代编程语言中非常重要的一部分，主机端将存储组织成以下结构：
一个程序不是每时每刻都需要访问其内存的任意部分，而是表现出局部性：如果一个地址被访问一次，那么在近期内它很可能还会被访问；如果一个地址被访问一次，那么它附近的地址也可能会在接下来被访问。比如指令的执行，在遇到跳转指令之前，指令总是连续存放的，而跳转指令有相当一部分是循环，循环中的指令都将被执行许多次。
基于局部性原理，现代CPU的储存体系被组织成多级结构，越靠近CPU的越快，每bit的造价也越高，因此容量也越低。一般而言，除磁盘以外，其它各级存储都不会保存所有数据，如内存使用请求调页机制将需要访问的页面调入内存，而缓存将正在访问的地址所在的行调入缓存。从这个过程也可以看出，对存储器的访问有一个最小单位，它不一定是机器的字长，比如内存页面的大小为4KB，则一次请求调页必然从磁盘调入至少4KB的数据，而如果缓存行的大小为64字节，则一次访存必然从内存读入64字节的数据。因此在设计程序时应当仔细考虑数据访问的顺序，否则缓存的频繁失效将显著降低程序的性能。
比如矩阵乘法：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
|
// C = A x B。dim：MxN = MxK * KxN
// C[i][j]为A的第i行乘B的第j列
for (int i = 0; i < M; ++i) {
for (int j = 0; j < N; ++j) {
C[i][j] = 0;
for (int l = 0; l < K; ++l)
C[i][j] += A[i][l] * B[l][j]; // B的访问缓存不友好
}
}
// C的第i行是以A的第i行各个元素为系数对B的各行进行线性组合的结果
for (int i = 0; i < M; ++i) {
for (int j = 0; j < N; ++j)
C[i][j] = 0;
for (int l = 0; l < K; ++l) {
for (int j = 0; j < N; ++j)
C[i][j] += A[i][l] * B[l][j]; // A、B、C都按行访问
}
}
|
CUDA的储存体系为多线程执行而设计，在储存体系的组织上与主机端有很大的不同，大致的结构为：
CUDA内存模型中的存储器类型比CPU丰富得多，有：寄存器、共享内存、本地内存、常量内存、纹理内存、全局内存、一级和二级缓存。其中一级缓存和二级缓存是不可编程的，其它存储器是程序员可以操作的。不同的储存器有不同的特点，使用场景也不同。
寄存器是GPU上速度最快的存储器，核函数中声明的一个没有其他修饰符的变量，通常存储在寄存器中。在核函数中声明的数组，如果访问时的索引都是编译期可以确定的常数，那么该数组也存储在寄存器中。比如：
1
2
3
4
5
6
7
|
__global__ void kernel() {
unsigned int id = threadIdx.x + blockIdx.x * blockDim.x;
int A[3];
A[0] = 1;
A[1] = 2;
A[2] = 3;
}
|
一个流式多处理器拥有数万个寄存器，但分配到每个线程后仍然不是太多，比如pascal架构的一个SM有65536个32位寄存器，但分配给2048个线程后每个线程只有32个寄存器，如果使用了double等长度为64位的变量则会占用两个寄存器。如果一个核函数需要的寄存器数量超过了硬件限制，就会使用本地内存代替多占用的寄存器，这将给性能带来非常不利的影响。
本地内存是每个线程独有的，用来存放：1、在编译时未知索引的数组，2、可能占用大量寄存器空间的结构体或者数组，3、超过硬件数量限制的变量(register spilling)。本地内存在逻辑上属于每个线程，但在存储体上无独立硬件，它和全局内存是同一个存储空间，因此速度慢、延迟高，使用时访存的优化和全局内存类似。
共享内存位于SM内部，在Fermi和Kepler架构中，共享内存和一级缓存是同一个硬件，默认配置下，一级缓存为16KB，共享内存为48KB，可以通过cudaDeviceSetCacheConfig等函数来设置它们的大小。从Pascal架构开始，共享内存已经独立出来，为64KB，前述的配置函数调用时不会出错，但也已经没有效果了。共享内存使用时通过__shared__来声明，它在线程块内共享，如：
1
2
3
|
__global__ void kernel() {
__shared__ double cache[1024]; // 同一线程块内的cache是同一块内存
}
|
共享内存在核函数的范围内声明，其生命周期和线程块的生命周期相同。共享内存是线程之间相互通信的基本方式，一个块内的线程通过共享内存中的数据可以相互合作。
全局内存是容量最大的内存，和主机端的内存对应，访存的模型也类似。全局内存通过函数cudaMalloc分配，通过cudaFree释放。一般在计算的时候将全局内存的指针传递给核函数，在全局内存上分配的空间存在于程序的整个生命周期，可被核函数的所有线程访问。需要注意，线程的执行不能跨线程块同步，因此多个线程修改全局内存的同一个位置可能出现问题，需要使用原子操作进行。和内存类似，访问全局内存时会将附近的数据缓存到二级缓存和一级缓存，缓存行的大小一般为128字节，也就是说一次内存事务将读入128字节的数据。需要注意，一个SM有16KB的一级缓存，最多要为2048个线程服务，缓存失效可能会非常快，读入的数据必须即使使用，否则必须重新读取。
常量内存是全局的，是一种专用内存，它的访问不占用全局内存的带宽和L1、L2缓存，其缓存位于SM内部，每个SM有64KB的常量缓存。常量内存的访问特点是硬件会将取到的数广播给线程束内一半的线程，它对于核函数来说是只读的，但是可以通过主机端对它进行更改。常量内存通过__constant__来声明，其声明应当放在全局作用域，对编译单元内的核函数可见，在主机端通过函数cudaMemcpyToSymbol来拷贝数据到常量内存。需要注意常量内存的访问要尽量满足线程束内多个线程访问同一个地址，因为以不同的地址访问常量内存将以串行的方式工作，性能可能还不如全局内存。
纹理内存的模型与常量内存类似，也是只读的，并在SM内部缓存。纹理内存的特点是多维的局部性，比如图片的访问一般不是按行来访问的，在访问了一个地址之后，接下来的访问位置很可能在图片上对应点的周围，因此纹理内存的缓存不是按行来缓存的，而是根据空间的局部性进行缓存。
CUDA的各种存储体的特性总结见下表：
全局内存和主机端的内存非常像，所有对全局内存的访问都会经过二级缓存，也可能会经过一级缓存，这取决于GPU的架构，比如计算能力3.x的设置默认不将数据缓存到一级缓存，但计算能力3.5和3.7的一些设备可以通过编译参数-Xptxas -dlcm=ca来启用用一级缓存。这些细节上的不同可以在Programming Guide的I. Compute Capabilities一节中查到。如果只使用二级缓存，一次内存访问将由一个32字节的内存事务来实现，如果使用了一级缓存，则一次内存访问将由一个128字节的内存事务来实现。
一行一级缓存是128个字节，如果线程束中的每个线程都请求一个4字节的值，总共就需要128字节的数据。刚好和一级缓存的一行大小相同，也是一次内存事务读入的字节数量。因此在设计优化程序时，需要考虑两个问题：1、对齐的内存访问；2、合并的内存访问。
当设备内存事务的第一个地址是用于事务服务的缓存粒度的偶数倍时(32字节的二级缓存或128字节的一级缓存)，就会出现对齐内存访问。运行非对齐的加载会造成带宽浪费。当一个线程束中全部的32个线程访问一个连续的内存块时，就会出现合并内存访问。对齐合并内存访问的理想状态是线程束从对齐内存地址开始访问一个连续的内存块。如下面两图所示，前者为对齐的内存访问，仅需要一次128字节的内存事务，而后者可能需要三次，或者前者需要4次32字节的内存事务而后者需要12次。
在最坏的情况下，线程束中线程访问分散于全局内存中的32个4字节的数据，完成此次访问需要32次内存事务，总共有31次内存事务的数据量未被使用。
GPU的缓存和CPU的缓存的作用有些许微妙的不同，CPU的缓存根据时间和空间的局部性优化访存，而GPU的二级缓存作用类似，但一级缓存几乎只能优化空间局部性而无法优化时间局部性，这是因为一个SM只有16KB的一级缓存，而最多可以运行2048个线程，访存的数据总量非常高，频繁访问一个在一级缓存的内存位置不会增加数据留在缓存中的概率。
如果访存确实无法对齐，或者无法合并，则更小的内存事务可以带来更高的访存效率，也就是说，不启用一级缓存反而比启用一级缓存效率更高！
对于数据的存储，一般有两种组织方案：结构体数组和数组结构体，比如二维平面上的点集：
1
2
3
4
5
6
7
8
9
10
11
12
13
|
// 结构体数组
struct Point {
double x;
double y;
};
Point points[1024];
// 数组结构体
struct Points {
double x[1024];
double y[1024];
};
Points points;
|
一般而言，大家倾向于第一种方案，它更符合我们对类型的认识，但它不利于访存，在一次内存事务中，可能只有一半的数据(x坐标或者y坐标)是有用的，而另一半的数据如果能在缓存失效前访问问题也不大，如果不能就会造成内存带宽的浪费。因此最好还是按照最第二个方案来设计程序。
在CUDA编程中，拆分任务是非常重要的一部分，但不能像CPU并行那样仅考虑任务计算负载的均衡，更重要的是要考虑到访存的性能，因为GPU的计算峰值非常高，计算往往不是性能瓶颈，而访存才是，这是正是GPU拥有类型如此丰富的存储器的原因。比如对于矩阵和向量的乘法，一个比较简单的想法是一个线程计算一个分量，一般在CPU上的并行就是这样做的：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
|
constexpr int M = 10000;
constexpr int N = 20000;
std::vector<double> mat(M * N);
std::vector<double> vec(N);
std::vector<double> result(M);
// generate mat and vec
tbb::parallel_for(tbb::blocked_range<int>(0, M),
[&mat, &vec, &result, N](const tbb::blocked_range<int> &r) -> void {
for (int i = r.begin(), i != r.end(); ++i) {
double t = 0;
for (int j = 0; j < N; ++j)
t += mat[i * N + j] * vec[j];
result[i] = t;
}
});
|
对于CPU而言没有问题，对内存的访问也是连续的地址，缓存的命中率非常高。如果在CUDA中也用一个线程计算一行，代码大致为：
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
|
constexpr int M = 10000;
constexpr int N = 20000;
__device__ inline double getMatElement(const double *mat, int row, int col, int N) {
return mat[row * N + col];
}
__global__ void MatMulVecKernel(
double *result,
const double *mat,
const double *vec,
const int M,
const int N
) {
unsigned int i = threadIdx.x + blockIdx.x * blockDim.x;
if (i < M) {
double t = 0;
for (int j = 0; j < N; ++j)
t += getMatElement(mat, i, j, N) * vec[j];
result[i] = 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 + block.x - 1) / block.x};
MatMulVecKernel<<<grid, block>>>(result, mat, vec, M, N);
// use result
cudaFree(mat);
cudaFree(vec);
cudaFree(result);
}
|
显然，在GPU的内存模型下，这样的计算方式性能非常低，因为同一个线程束内的32个线程访问矩阵相邻的32行，每次访存的地址相差N * sizeof(double)，在N大于等于16时需要32次内存事务才能取得所有数据，这是前面所说的最坏情况。为了合并访存，应当至少使用32个线程来处理矩阵的一行，且相邻32个线程处理相邻的数据，比如使用32个线程计算一个分量：
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);
}
|
两者在计算能力为3.5为k20m上的时间对比为：
uncoalesced: 0.105597 second
coalesced: 0.0119435 second
性能差距接近10倍。
完整的测试代码：
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
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
|
#include <iostream>
#include <iomanip>
#include <vector>
#include <random>
#include <chrono>
#include <cuda_runtime.h>
#include <curand.h>
#define CHECK(call) \
do { \
const cudaError_t error = call; \
if (error != cudaSuccess) \
{ \
fprintf(stderr, "Error: %s:%d, ", __FILE__, __LINE__); \
fprintf(stderr, "code: %d, reason: %s\n", error, \
cudaGetErrorString(error)); \
exit(1); \
} \
} while (0)
__global__ void matrixVectorMulKernel0(
double *out_vec,
const double *mat,
const double *in_vec,
const int M,
const int N
) {
unsigned int row = threadIdx.x + blockIdx.x * blockDim.x;
if (row < M) {
double t = 0;
for (int col = 0; col < N; ++col)
t += mat[row * N + col] * in_vec[col];
out_vec[row] = t;
}
}
__global__ void matrixVectorMulKernel(
double *out_vec,
const double *mat,
const double *in_vec,
volatile double *cache,
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[threadIdx.x];
}
}
void generateMatrix(double *mat, const int M, const int N) {
curandGenerator_t generator;
CHECK_CURAND(curandCreateGenerator(&generator, CURAND_RNG_PSEUDO_DEFAULT));
CHECK_CURAND(curandSetPseudoRandomGeneratorSeed(generator, std::random_device{}()));
CHECK_CURAND(curandGenerateUniformDouble(generator, mat, M * N));
}
void generateVector(double *vec, const int N) {
curandGenerator_t generator;
CHECK_CURAND(curandCreateGenerator(&generator, CURAND_RNG_PSEUDO_DEFAULT));
CHECK_CURAND(curandSetPseudoRandomGeneratorSeed(generator, std::random_device{}()));
CHECK_CURAND(curandGenerateUniformDouble(generator, vec, N));
}
int main() {
std::chrono::steady_clock::time_point start, end;
std::chrono::duration<double, std::ratio<1>> elapsed;
const int M = 10000;
const int N = 20000;
double *mat = nullptr; // M x N
double *in_vec = nullptr; // N x 1
double *out_vec = nullptr; // M x 1
CHECK(cudaMalloc(&mat, M * N * sizeof(double)));
CHECK(cudaMalloc(&in_vec, N * sizeof(double)));
CHECK(cudaMalloc(&out_vec, M * sizeof(double)));
// 生成矩阵和向量
start = std::chrono::steady_clock::now();
generateMatrix(mat, M, N);
generateVector(in_vec, N);
CHECK(cudaDeviceSynchronize());
end = std::chrono::steady_clock::now();
elapsed = end - start;
std::cout << "generate matrix and vector: " << elapsed.count() << " second" << std::endl;
dim3 block;
dim3 grid;
// 不合并的访存
block.x = 512;
grid.x = (M + block.x - 1) / block.x;
start = std::chrono::steady_clock::now();
matrixVectorMulKernel0<<<grid, block>>>(out_vec, mat, in_vec, M, N);
CHECK(cudaDeviceSynchronize());
end = std::chrono::steady_clock::now();
elapsed = end - start;
std::cout << "uncoalesced: " << elapsed.count() << " second" << std::endl;
// 合并的访存
block.x = 512;
grid.x = (M * 32 + block.x - 1) / block.x;
double *cache = nullptr;
CHECK(cudaMalloc(&cache, block.x * grid.x * sizeof(double)));
start = std::chrono::steady_clock::now();
matrixVectorMulKernel<<<grid, block>>>(out_vec, mat, in_vec, cache, M, N);
CHECK(cudaDeviceSynchronize());
end = std::chrono::steady_clock::now();
elapsed = end - start;
std::cout << "coalesced: " << elapsed.count() << " second" << std::endl;
cudaFree(mat);
cudaFree(in_vec);
cudaFree(out_vec);
cudaFree(cache);
return 0;
}
|