CUDA 矩阵乘法优化

矩阵乘法
　　为了单纯起见，我们这里以方形的矩阵为例子。基本上，假设有两个矩阵 A 和 B，则计算 AB = C 的方法如下：

for(j
= 0; j< n; j++)
{

C[i][j] =0;

for(k=
0; k< n; k++) {

C[i][j] += A[i][k]* B[k][j];

}

}

}

一开始，我们先准备好产生数据、设定 CUDA 等等的工作：

int main()

{

float*a,
*b,
*c,*d;

int n=
1000;

if(!InitCUDA())return
0;

a = (float*) malloc(sizeof(float)*
n * n);

b = (float*) malloc(sizeof(float)*
n * n);

c = (float*) malloc(sizeof(float)*
n * n);

d = (float*) malloc(sizeof(float)*
n * n);

srand(0);

matgen(a, n, n);

matgen(b, n, n);

clock_t time = matmultCUDA(a, n, b, n, c, n, n);

matmult(a, n, b, n, d, n, n);

compare_mat(c, n, d, n, n);

double sec= (double) time/
CLOCKS_PER_SEC;

printf("Time used: %.2f (%.2lf GFLOPS)\n", sec,

2.0* n
* n
* n/ (sec
* 1E9));

return0;

}

InitCUDA 函式和第一个 CUDA 程序一样，可以直接参考前面的文章。以下是上面用到的一些其它的函式：

　　产生矩阵：

void matgen(float* a,int
lda, int n)

{

int i, j;

for(i=
0; i< n; i++) {

for(j=
0; j< n; j++) {

a[i * lda+ j]
= (float) rand()/ RAND_MAX
+

(float) rand()/ (RAND_MAX*
RAND_MAX);

}

}

}

这个函式只是利用随机数生成器把矩阵填满 0 ~ 1 之间的数字。特别注意到因为 C 语言中无法声明变动大小的二维矩阵，所以我们使用 i * lda + j 的方式。

　　进行矩阵乘法：

void matmult(constfloat*
a,int lda,
const
float* b,int ldb,

float* c,int ldc,
int n)

{

int i, j, k;

for(i=
0; i< n; i++) {

for(j=
0; j< n; j++) {

double t=
0;

for(k=
0; k< n; k++) {

t += a[i* lda
+ k]
* b[k* ldb
+ j];

}

c[i * ldc+ j]
= t;

}

}

}

这是以
CPU 进行矩阵乘法、用来进行验证答案正确与否的程序。特别注意到它用 double 来储存暂时的计算结果，以提高精确度。

　　验证结果：

void compare_mat(constfloat*
a,int lda,

constfloat* b,int
ldb, int n)

{

float max_err=
0;

float average_err=
0;

int i, j;

for(i=
0; i< n; i++) {

for(j=
0; j< n; j++) {

if(b[i* ldb
+ j]
!=0) {

float err= fabs((a[i* lda
+ j]
-

b[i * ldb+ j])
/ b[i
* ldb+ j]);

if(max_err< err) max_err= err;

average_err += err;

}

}

}

printf("Max error: %g Average error: %g\n",

max_err, average_err / (n* n));

}

这个函式计算两个矩阵的最大相对误差和平均相对误差，并把结果印出来。

　　最后是 CUDA 的矩阵乘法的部份：

#define NUM_THREADS 256

clock_t matmultCUDA(constfloat* a,int
lda,

constfloat* b,int
ldb, float* c,int ldc,
int n)

{

float*ac,
*bc,
*cc;

clock_t start, end;

start = clock();

cudaMalloc((void**)&ac,
sizeof(float)* n
* n);

cudaMalloc((void**)&bc,
sizeof(float)* n
* n);

cudaMalloc((void**)&cc,
sizeof(float)* n
* n);

cudaMemcpy2D(ac, sizeof(float)* n, a,
sizeof(float)* lda,

sizeof(float)* n, n, cudaMemcpyHostToDevice);

cudaMemcpy2D(bc, sizeof(float)* n, b,
sizeof(float)* ldb,

sizeof(float)* n, n, cudaMemcpyHostToDevice);

int blocks= (n
+ NUM_THREADS
- 1)/ NUM_THREADS;

matMultCUDA<<<blocks* n, NUM_THREADS>>>

(ac, n, bc, n, cc, n, n);

cudaMemcpy2D(c, sizeof(float)* ldc, cc,
sizeof(float)* n,

sizeof(float)* n, n, cudaMemcpyDeviceToHost);

cudaFree(ac);

cudaFree(bc);

cudaFree(cc);

end = clock();

return end- start;

}

这个函式相当单纯，就是在显卡内存中配置存放矩阵的内存，然后把主内存中的矩阵数据复制到显卡内存上。不过，因为我们的矩阵乘法函式可以指定
pitch(即 lda、ldb、和 ldc)，所以如果用一般的 cudaMemcpy 函式来复制内存的话，会需要每个 row 都分开复制，那会需要呼叫很多次 cudaMemcpy 函式，会使效率变得很差。因此，在这里我们用了一个新的 cudaMemcpy2D 函式，它是用来复制二维数组，可以指定数组的 pitch。这样就可以透过一次函数调用就可以了。

　　进行计算的 kernel 如下：

__global__ staticvoid matMultCUDA(constfloat*
a, size_t lda,

constfloat* b, size_t ldb,float*
c, size_t ldc,int n)

{

constint tid
= threadIdx.x;

constint bid
= blockIdx.x;

constint idx
= bid
* blockDim.x+ tid;

constint row
= idx
/ n;

constint column
= idx
% n;

int i;

if(row< n
&& column
< n) {

float t=
0;

for(i=
0; i< n; i++) {

t += a[row* lda
+ i]
* b[i* ldb
+ column];

}

c[row * ldc+ column]
= t;

}

}

这个函式一开始先从 bid 和 tid 计算出这个 thread 应该计算的 row 和 column，在判断 row 和 column 在范围内之后，就直接进行计算，并把结果写到 c 矩阵中，是非常单纯的函式。

　　在 GeForce 8800GT 上实际执行的结果如下：

　　Max error: 2.01484e-006 Average error: 3.36637e-007

　　Time used: 1.1560 (1.73 GFLOPS)

　　可以看到两个问题：

　　很明显的，执行效率相当低落。

　　最大相对误差偏高。理想上应该要低于 1e-6。

　　计算结果的误差偏高的原因是，在
CPU 上进行计算时，我们使用 double(即 64 bits 浮点数)来累进计算过程，而在 GPU 上则只能用 float(32 bits 浮点数)。在累加大量数字的时候，由于累加结果很快会变大，因此后面的数字很容易被舍去过多的位数。

　　由于 CUDA 的浮点数运算，在进行加、减、乘法时是符合 IEEE 754 规定的精确度的，因此，我们可以利用 Kahan's Summation Formula 来提高精确度。把程序改成：

if(row < n&& column
< n) {

float t=
0;

float y=
0;

for(i=
0; i< n; i++) {

float r;

y -= a[row* lda
+ i]
* b[i* ldb
+ column];

r = t- y;

y = (r- t)
+ y;

t = r;

}

}

修改后的程序的执行结果是：

　　Max error: 1.19209e-007 Average error: 4.22751e-008

　　Time used: 1.1560 (1.73 GFLOPS)

　　可以看到相对误差有很大的改善，效率则没什么变化。

　　由于 Kahan's Summation Formula 需要的运算量提高，但是效率却没有什么改变，可以看出这个 kernel 主要的瓶颈应该是在内存的存取动作上。这是因为有大量的内存读取是重复的。例如，矩阵 a 的一个 row 在每次进行计算时都被重复读入，但这是相当浪费的。这样的计算方式，总共需要读取 2*n3 次内存。如果让一个 row 只需要读入一次的话，就可以减到为 n3+n2 次。

　　第一个改良

　　和我们的第一个 CUDA 程序一样，我们可以利用 shared memory 来储存每个 row 的数据。不过，因为只有同一个 block 的 threads 可以共享 shared memory，因此现在一个 row 只能由同一个 block 的 threads 来进行计算。另外我们也需要能存放一整个 row 的 shared memory。因此，把先把呼叫 kernel 的部份改成：

　　matMultCUDA<<>>

　　(ac, n, bc, n, cc, n, n);

　　kernel 的部份则改成：

__global__ staticvoid matMultCUDA(constfloat*
a, size_t lda,

constfloat* b, size_t ldb,float*
c, size_t ldc,int n)

{

extern __shared__float data[];

constint tid
= threadIdx.x;

constint row
= blockIdx.x;

int i, j;

for(i= tid; i
< n; i
+= blockDim.x) {

data[i] = a[row* lda
+ i];

}

__syncthreads();

for(j= tid; j
< n; j
+= blockDim.x) {

float t=
0;

float y=
0;

for(i=
0; i< n; i++) {

float r;

y -= data[i]* b[i
* ldb
+ j];

r = t- y;

y = (r- t)
+ y;

t = r;

}

c[row * ldc+ j]
= t;

}

}

第一个部份先把整个 row 读到 shared memory 中，而第二个部份则进行计算，并没有太大的变化。主要的差别是现在一个 row 只由一个 block 进行计算。

　　在 GeForce 8800GT 上，执行的结果是：

　　Max error: 1.19209e-007 Average error: 4.22751e-008

　　Time used: 0.4220 (4.74 GFLOPS)

　　很明显的，计算的结果并没有改变，不过速度则提高了超过一倍。虽然如此，但是这样的效率仍不尽理想，因为理论上 GeForce 8800GT 有超过 300GFLOPS 的运算性能。即使是把 Kahan's Summation Formula 所需要的额外运算考虑进去，这样的效率仍然连理论最大值的十分之一都不到。

　　会有这样的结果，原因其实还是同样的：对内存的存取次数太多了。虽然现在 A 矩阵的 row 的数据已经不再需要重复读取，但是 B 矩阵的 column 的数据仍然一直被重复读取。

　　另一个问题比较不是那么明显：对 B 矩阵的读取，虽然看起来不连续，但实际上它是连续的。这是因为不同的 thread 会读取不同的 column，因此同时间每个 thread 读取的各个 column 加起来，就是一个连续的内存区块。那么，为什么效率还是不佳呢?这是因为，GPU 上的内存控制器，从某个固定的倍数地址开始读取，才会有最高的效率(例如 16 bytes 的倍数)。由于矩阵大小并不是 16 的倍数(这里使用的是 1000x1000 的矩阵)，所以造成效率不佳的情形。

　　要解决这个问题，我们可以在 cudaMalloc 的时候稍微修改一下，让宽度变成 适当的倍数就可以了。但是，适当的倍数是多少呢?幸运的是，我们并不需要知道这些细节。CUDA 提供了一个 cudaMallocPitch 的函式，可以自动以最佳的倍数来配置内存。因此，我们可以把 cudaMalloc 的部份改成：

size_t pitch_a, pitch_b, pitch_c;

cudaMallocPitch((void**)&ac,
&pitch_a,
sizeof(float)* n, n);

cudaMallocPitch((void**)&bc,
&pitch_b,
sizeof(float)* n, n);

cudaMallocPitch((void**)&cc,
&pitch_c,
sizeof(float)* n, n);

cudaMallocPitch 函式会以适当的倍数配置内存，并把配置的宽度传回。因此，在把矩阵复制到显卡内存上时，要使用它传回的宽度：

cudaMemcpy2D(ac, pitch_a, a, sizeof(float)
* lda,

sizeof(float)
* n, n, cudaMemcpyHostToDevice);

cudaMemcpy2D(bc, pitch_b, b, sizeof(float)
* ldb,

sizeof(float)
* n, n, cudaMemcpyHostToDevice);

呼叫 kernel 的部份也需要修改：

　　matMultCUDA<<>>

　　(ac, pitch_a / sizeof(float), bc, pitch_b / sizeof(float),

　　cc, pitch_c / sizeof(float), n);

　　同样的，把计算结果复制回到主内存时，也要使用传回的宽度值：

　　cudaMemcpy2D(c, sizeof(float) * ldc, cc, pitch_c,

　　sizeof(float) * n, n, cudaMemcpyDeviceToHost);

　　这样就完成了。Kernel 部份则不需要修改。

　　这样的修改有多大的效果呢?在 GeForce 8800GT 上执行，结果如下：

　　Max error: 1.19209e-007 Average error: 4.22751e-008

　　Time used: 0.1250 (16.00 GFLOPS)

　　可以看到，执行速度又再大幅提高了三倍多!而这只是把内存的配置方式稍微修改一下而已。

　　虽然执行速度提高了很多，但是，和前面提到的理论值相比，其实还是有相当的差距。这是因为，前面也提到过，这样的做法需要 n3+n2 次的内存读取，和 n2 次的内存写入动作。由于 n = 1000，每个数字的大小是 32 bits，所以总共的内存存取数据量约为 4GB。除以实际执行的时间 0.125 秒，得到的带宽数值是约 32GB/s，这已经接近 GeForce 8800GT 显卡内存的带宽了。由于我们计算时间的时候，把配置内存、以及数据的复制动作也计算进去，因此实际上花费在 kernel 的时间是更短的(约
0.09 秒)。因此，可以很明显的看出，这个程序的效率，是受限于内存带宽的。

　　进一步的改良

　　上一节的结论显示出，矩阵乘法的程序，效率是受限于内存带宽的。那有没有办法降低内存的存取次数呢?答案当然是有的，不然就不会有这一节了 :)

　　要进一步降低内存带宽的使用，可以注意到，在上一节的方法中，虽然 A 矩阵的存取次数被减至最低，但是 B 矩阵的存取次数并没有减少。这是因为我们只将 A 矩阵的 row 加载到 shared memory 中，但是 B 矩阵的 column 也是有被重复使用的。理想上应该也可以避免重复加载才对。不过，由于 B 矩阵的 column 使用的时机，和 A 矩阵的 row 是不同的，所以并不能直接这样做。

　　解决方法是 "blocking"。也就是把整个矩阵乘法的动作，切割成很多小矩阵的乘法。例如，要计算 C 矩阵的 (0, 0) ~ (15, 15) 的值，可以把它想成：

　　A(0~15, 0~15) * B(0~15, 0~15) + A(0~15,16~31) * B(16~31, 0~15)

　　+ A(0~15, 32~47) * B(32~47, 0~15) + ...

　　这样一来，我们就可以把两个小矩阵加载到 shared memory，则小矩阵本身的乘法就不需要再存取任何外部的内存了!这样一来，假设小矩阵的大小是 k，则实际上需要的内存存取次数就会变成约 2k2(n/k)3 = 2n3/k。

　　由于目前 CUDA 每个 block 的 thread 数目最多是 512，因此 k = 16 似乎是一个相当理想的数字(共 256 个 threads)。因此，对于一个 n = 1000 的矩阵来说，我们可以把内存存取的量减少到约 500MB，也就是上一节的存取量的 1/8。理论上，这样应该可以让效率提高八倍(假设没有遇到别的瓶颈)。

　　为了方便进行区块的计算，我们让每个 block 有 16x16 个 threads，再建立 (n/16)x(n/16) 个 blocks。把呼叫 kernel 的地方改成：

int bx
= (n + BLOCK_SIZE
- 1)
/ BLOCK_SIZE;

dim3 blocks(bx, bx);

dim3 threads(BLOCK_SIZE, BLOCK_SIZE);

matMultCUDA<<<blocks, threads>>>(ac, pitch_a
/ sizeof(float),

bc, pitch_b /
sizeof(float), cc, pitch_c
/ sizeof(float), n);

　　BLOCK_SIZE 则是定义成 16。dim3 是 CUDA 的一种数据型态，表示一个 3D 的向量。在这里，我们透过 dim3 来建立 16x16 个 threads 的 block，和 (n/16)x(n/16) 个 blocks。

　　Kernel 程序的部份，则改成：

__global__ static
void matMultCUDA(const
float* a, size_t lda,

const
float* b, size_t ldb,
float* c, size_t ldc,
int n)

{

__shared__ float matA[BLOCK_SIZE][BLOCK_SIZE];

__shared__ float matB[BLOCK_SIZE][BLOCK_SIZE];

const
int tidc
= threadIdx.x;

const
int tidr
= threadIdx.y;

const
int bidc
= blockIdx.x
* BLOCK_SIZE;

const
int bidr
= blockIdx.y
* BLOCK_SIZE;

int i, j;

float results
= 0;

float comp
= 0;

for(j
= 0; j
< n; j
+= BLOCK_SIZE) {

if(tidr
+ bidr
< n && tidc
+ j
< n) {

matA[tidr][tidc] = a[(tidr
+ bidr)
* lda + tidc
+ j];

}

else {

matA[tidr][tidc] =
0;

}

if(tidr
+ j
< n && tidc
+ bidc
< n) {

matB[tidr][tidc] = b[(tidr
+ j)
* ldb + tidc
+ bidc];

}

else {

matB[tidr][tidc] =
0;

}

__syncthreads();

for(i
= 0; i
< BLOCK_SIZE; i++) {

float t;

comp -= matA[tidr][i]
* matB[i][tidc];

t = results
- comp;

comp = (t
- results)
+ comp;

results = t;

}

__syncthreads();

}

if(tidr
+ bidr
< n && tidc
+ bidc
< n) {

c[(tidr + bidr)
* ldc
+ tidc + bidc]
= results;

}

}

注意到因为我们现在使用 16x16 的 threads，因此 threadIdx 变量可以取得 threadIdx.x 和 threadIdx.y，范围分别是 0 ~ 15。blockIdx.x 和 blockIdx.y 变量也是同样的情形，范围分别是 0 ~ n/16。

　　在程序中，因为矩阵的大小不一定会是 16 的倍数，因此需要使用 if 判断式检查是否超出矩阵范围。

　　这个版本在 GeForce 8800GT 上的执行结果如下：

　　Max error: 1.19209e-007 Average error: 4.22751e-008

　　Time used: 0.0780 (25.64 GFLOPS)

　　速度虽然提高了，但是似乎并没有达到预期中的八倍。当然，前面提到过，我们在计算时间时，把一些复制内存、配置内存的动作也计算在内，这些动作的时间并不会缩短。实际上 kernel 的运行时间，大约是 0.053 秒左右(约略相当于 38GFLOPS)，比上一节的版本快了将近一倍。

　　如果这一版程序已经不再限于内存带宽，那为什么没有达到预期的效率呢?这是因为这一版程序已经是限于指令周期了。除了使用 Kahan's Summation Formula 会需要更多的运算之外，程序中也有大量计算矩阵地址的乘法等等，这都会需要花费运算资源。另外，那些用来判断超出矩阵范围的 if 判断式，也会有一定的影响。

　　要把那些 if 判断式去掉，有一个方法是，在配置内存时，就配置成 16 的倍数，并在复制矩阵到显卡内存之前，先将它清为 0。如下所示：

int newn = ((n
+ BLOCK_SIZE
- 1)
/ BLOCK_SIZE)
* BLOCK_SIZE;

cudaMallocPitch((void**)
&ac,
&pitch_a,

sizeof(float)
* newn, newn);

cudaMallocPitch((void**)
&bc,
&pitch_b,

sizeof(float)
* newn, newn);

cudaMallocPitch((void**)
&cc,
&pitch_c,

sizeof(float)
* newn, newn);

cudaMemset(ac, 0, pitch_a
* newn);

cudaMemset(bc, 0, pitch_b
* newn);
　　这样一来，我们就可以把 kernel 中的 if 判断式都移除了：

__global__ static
void matMultCUDA(const
float* a, size_t lda,

const
float* b, size_t ldb,
float* c, size_t ldc,
int n)

{

__shared__ float matA[BLOCK_SIZE][BLOCK_SIZE];

__shared__ float matB[BLOCK_SIZE][BLOCK_SIZE];

const
int tidc
= threadIdx.x;

const
int tidr
= threadIdx.y;

const
int bidc
= blockIdx.x
* BLOCK_SIZE;

const
int bidr
= blockIdx.y
* BLOCK_SIZE;

int i, j;

float results
= 0;

float comp
= 0;

for(j
= 0; j
< n; j
+= BLOCK_SIZE) {

matA[tidr][tidc] = a[(tidr
+ bidr)
* lda + tidc
+ j];

matB[tidr][tidc] = b[(tidr
+ j)
* ldb + tidc
+ bidc];

__syncthreads();

for(i
= 0; i
< BLOCK_SIZE; i++) {

float t;

comp -= matA[tidr][i]
* matB[i][tidc];

t = results
- comp;

comp = (t
- results)
+ comp;

results = t;

}

__syncthreads();

}

c[(tidr + bidr)
* ldc
+ tidc + bidc]
= results;

}
　　这个版本的执行结果是：

　　Max error: 1.19209e-007 Average error: 4.22751e-008

　　Time used: 0.0780 (25.64 GFLOPS)

　　似乎没有改善。不过，实际上 kernel 的运行时间已经减少到 0.042 秒(约略相当于 48GFLOPS)。

　　结论

　　有些读者可能会想，如果把 block 再变得更大(例如 32x32)是否会有帮助呢?当然，由于最后的程序已经不再是受限于内存带宽(在 0.042 秒内存取 500MB 的数据约相当于 12GB/s 的带宽)，所以把 block 再加大并不会有帮助了。而且，由于一个 block 内的 thread 数目最多只能到 512 个，将 block 变大也会造成很多额外负担。而且 shared memory 的大小也有限制(GeForce 8800GT 的 shared memory 大小限制是 16384
bytes)，所以也不能任意增加 block 的大小。