所有产品

和记娱乐cuda入门:如何进行矩阵乘法优化

  【IT168技术】前面介绍的计算平方和的程序,似乎没有什么实用价值。所以我们的第二个 CUDA 程序,要做一个确实有(某些)实用价值的程序,也就是进行矩阵乘法。而且,这次我们会使用浮点数。

  虽然矩阵乘法有点老套,不过因为它相当简单,而且也可以用来介绍一些有关 CUDA 的有趣性质。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  虽然执行速度提高了很多,但是,和前面提到的理论值相比,其实还是有相当的差距。这是因为,前面也提到过,这样的做法需要 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) 的值,可以把它想成:

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

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

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

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

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

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

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