缓存结构GPU矩阵乘法算法的自动优化

李晓雯　崔翔　殷瑞杰　刘强

摘 要： 讨论在Fermi结构GPU使用CUDA对GEMM（单精度和双精度）算法进行优化，以及Fermi体系结构的新特性（如缓存）对性能的影响。GPU缓存一方面可以提高处理器在运行时数据访问的局部性，另一方面使得代码性能对与性能相关算法参数的依赖变得不可预测。自动优化技术可以用来解决这一问题。自动优化的SGEMM和DGEMM代码在Tesla C2050 GPU上达到了563GFlops和253GFlops的性能。代码使用CUDA和C语言进行实现，未进行二进制代码级别的优化。

关键词： GPU程序设计； 矩阵乘法； 自动优化； GEMM模板

中图分类号： TN40； TP312 文献标识码： A 文章编号： 1004?373X（2014）10?0137?04

Abstract： Automatic optimization of GEMM algorithm on Fermi GPU and the impact of Fermi GPUs architectural features on performance are discussed in this paper. The cache on GPU can not only improve the data access locality of processors， but also make the code performance relying on relative algorithm parameters unpredictable. Auto?tuning can be used to solve this problem. Auto?tuned SGEMM and DGEMM codes achieve 563 GFlops and 253 GFlops respectively on C2050 GPU. The codes are implemented by using CUDA and C language， but the optimization on the binary level is not involved.

Keywords： GPU programming； matrix multiplication； automatic optimization； GEMM template

0 引 言

Fermi是Nvidia公司支持CUDA编程模型的新一代GPU。与GT200体系结构相比，新型的Tesla 2050 GPU具有一些新的特性：如增强的双精度浮点性能、L1/L2 缓存结构、更多的寄存器、更大的共享存储器、ECC支持和更快的原子操作[1?3]。

由于Tesla 2050和GT200使用同样的编程模型，因此程序员期望在GT200上得到良好性能优化的代码也一样可以在Tesla 2050上取得良好的性能。实际上，程序员依然需要调整他们运行在GT200上的代码以在Tesla 2050上取得最高的性能。在Tesla 2050上，虽然每一个MP的寄存器文件被加倍了，但由于每一个MP的核数由8增加为32，因此实际上每个线程可用的寄存器数目实际上减少了一半。这使得程序员需要更加注意寄存器的使用效率。新增加的缓存结构一方面带来了运行时数据访问局部性的好处，另一方面也增加了代码性能的不可预测性。程序员依然需要了解GPU的硬件特性以得到高效的GPU代码。

自动调节技术是在复杂和不可预测的体系结构上得到近似最优代码的一种实用的技术。使用自动调节技术得到的SGEMM和DGEMM代码在Tesla 2050分别达到了563 GFlops和253 GFlops的速度，相对于CUBLAS 3.0分别具有1.7倍和1.6倍的加速[4]。

1 Fermi的新特性

1.1 L1/L2缓存

与GT200体系结构相比，Fermi增加了L1/L2缓存以提高访问设备存储器的性能，如图1所示。程序员对L1缓存的使用可以进行控制：64 KB的片上存储可以被用于L1缓存或共享存储器，而16 KB或48 KB存储分别用于L1缓存或共享存储器（抑或相反）可以在每次内核调用时进行控制。使用到局部存储器的内核代码可以从新增的L1缓存中受益。除了L1缓存，Fermi还提供768 KB的L2缓存。CUDA程序设计模型的原有特点是暴露硬件体系结构使得程序员对代码的性能可以进行良好的控制，而缓存的引入却使得CUDA代码的行为和性能变得难以预测。甚至CUDA编程手册[3]都建议程序员通过实验的方法来确定L1缓存或共享存储器的配置问题。对于L2缓存，一个简单的使得代码受益的方法是保证访问相同设备存储器地址空间的线程块被连续的调度；这可以通过将blockIdx.x和blockIdx.y变量进行对调实现。

考虑到缓存效果的不可预测性，自动优化技术可以用来得到近似最优的CUDA代码。首先将算法实现的代码进行参数模板化，通过选择不同的参数组合来自动得到具有良好性能的代码。

1.2 寄存器文件

在文献[4?5]中提到，相对于共享存储器而言，在设计算法时应该优先选择使用寄存器以得到良好的性能。在GT200体系结构上，使用寄存器间的MAD指令可以达到98%的理论运算性能。在Fermi发布之前，程序员都期望在新的GPU体系结构上，寄存器文件的大小能被增加，从而使得先前的代码在Fermi上能够取得更优的性能。

与GT200体系结构相比，Fermi上每个多处理器的寄存器文件的大小由16 KB增加为32 KB；而与此同时，每个多处理器的微核数目从8增加到32。这意味着在Fermi体系结构上，每个多处理器的微核的可用寄存器实际上是减少了。在GT200体系结构上，每个多处理器需要256个活动线程以掩盖指令流水线的延迟，而在Fermi体系结构上，则需要更多的活动线程来掩盖指令流水线的延迟。

1.3 32/64位设备代码

在Fermi体系结构上，如果代码按照64?bit的模式编译，则CUDA编译器会将CPU代码和设备代码都编译成为64?bit的目标代码。在这种情况下，设备代码中的指针变量将会占用多出一倍的寄存器空间。由于Tesla C2050的设备存储器容量不超过4 GB（在增加ECC的情况下可用的设备存储器空间只有2.625 GB），因此完全没有必要在设备代码中使用64 b的指针。因此，在本文实现的Fermi体系结构上的GEMM代码中，CPU代码和设备代码总是被分别编译的。

1.4 设备存储器访问

在GT200体系结构上，对设备存储器的访问是按照半个warp的单位来进行处理的，而在Fermi体系结构上，是按照一个warp的单位来进行处理的。因此，程序员需要调整内核调用时的维度设置。对于具有两个维度的线程块，其x维度的大小应该是warp大小的整倍数，而非半个warp大小的整倍数，从而使得每一个warp在访问设备存储器时可以得到较高的性能。

1.5 Bank冲突

在GT200体系结构上，共享存储器具有16个bank，而且对其的访问是按照半个warp的单位来进行处理的，而在Fermi体系结构上，共享存储器具有32个bank，而且对其的访问是按照整个warp的单位来进行处理的。每一个bank的大小为32位。

1.6 多内核并行执行

Fermi体系结构支持多个内核代码的并行执行，使得不同应用上下文的内核代码可以同时在一个GPU上运行；这样，多个小的内核代码可以共同利用一个GPU上的计算资源。这也是Fermi体系结构的新特性，但在本文自动优化的GEMM代码中并未使用到。

2 自动优化的GEMM代码

2.1 GEMM代码模板

文献[4]描述了在GTX280实现的达到393 Gflops性能的SGEMM内核代码，在此，依然使用该代码作为实现自动优化的代码模板。在该实现中，一个Csub被一个线程块进行计算，根据线程块中线程数目的多少，一个线程可以计算Csub的半列或多列元素。例如，对于m= 16和n=64，而线程块具有16×4个线程，则每一个线程将计算一整列Csub的元素；如果线程块具有16×8个线程，则每一个线程将计算半列Csub的元素。与文献[4]中的实现类似，每一个线程在使用线程ID计算出访问矩阵A、B和C的指针位置之后，进入一个循环。在每一轮循环中，一个线程块从设备存储器中读入一个Asub的数据到共享存储器中，之后，又一个内层的循环被执行：在每一轮内层循环中，一个线程从矩阵B中读入一个或多个元素，把这些元素与共享存储器中相应的数据做计算，将结果累加到Csub对应的寄存器中。最后，每个线程将其计算的Csub的数据写回到设备存储器中。

尽管与文献[4]中的GEMM实现具有类似的代码结构，但是为了满足缓存的友好性，对该代码模板进行了一个重要的修正。在此计算中，矩阵B总是位于设备存储器中，一列线程块总是需要读取一列Bsub的数据，如图2（a）所示。如果线程块能够按照列优先的顺序调度，则可以达到更好的缓存命中效果。因此，在此GEMM模板中，blockIdx.x和blockIdx.y变量的顺序被对调，从而达到如图2（b）所示的设备存储器访问效果。这样，在进行这个转置之后，线程块被以较优的缓存命中效果的方式调度。在后面给出的性能测试中，将对比这个对调进行或不进行的性能结果。

2.2 对代码模板进行自动优化

设计的自动调节程序根据代码模板生成代码并测试其性能结果。在代码模板中，5个参数（m，k，n，tx和ty）决定代码模板的执行行为；此外，引入2个额外的参数，一个用来决定线程块的维度是否对调，来测试缓存的效果，另一个确定L1缓存和共享存储器的配置比例。因此，整个代码模板的行为是由7个参数确定的。

图3显示经过优化得到的矩阵大小为2 048时的DGEMM代码的实例。注意：在该代码中，blockIdx.x和blockIdx.y变量通过C语言的宏进行了对调，而相应的调用代码也做了相应的改变。在此代码中，m=8，k= 64，n=1 024，tx=64，ty=8，因此一个线程负责计算Csub中的两列数据。

3 性能测试

4 结 语

为了在Fermi体系结构上书写高效的代码，程序员需要很好地了解Fermi体系结构的新特性，以及它们是如何影响程序的性能的。对于Fermi体系结构，程序员尤其要关心缓存对性能的影响。对于Fermi这种复杂的行为和性能难以预测的硬件体系结构，自动优化技术不失为得到高性能代码的一种实用技术。

参考文献

[1] NVIDIA Corp. Whitepaper： NVIDIA's next generation CUDA compute architecture [R/OL]. [2012?05?18]. http：//www. insidehpc.com.

[2] NVIDIA Corp. Tuning CUDA applications for Fermi [R/OL]. [2011?05?03]. http：// www. people.maths.ox.ac.uk/gilesm/cuda/doc/Fermi_Tuning_Guide.

[3] NVIDIA Corp. CUDA compute unified device architecture， programming guide， Version 3.0 [R/OL]. [2010?05?03]. http：// www. mohamedfahmed.wordpress.com

[4] CUI Xiang， CHEN Yi?feng， MEI Hong， et al. Auto?tuning GEMM for GPGPU with Cache [C]// Proceedings of 2010 IEEE 16th International Conference on Parallel and Distributed Systems （ICPADS）. Shanghai， China： IEEE， 2010： 237?242.

[5] VOLKOV V， DEMMEL J W. Benchmarking GPUs to tune dense linear algebra [C]// Proceedings of 2008. International Conference for High Performance Computing， Networking， Storage and Analysis Austin， TX： [s.n.]， 2008： 1?11.

[6] 李晓雯，崔翔.GPU矩阵乘法和FFT算法的性能优化[J].现代电子技术，2013，36（4）：80?84.

1.3 32/64位设备代码

在Fermi体系结构上，如果代码按照64?bit的模式编译，则CUDA编译器会将CPU代码和设备代码都编译成为64?bit的目标代码。在这种情况下，设备代码中的指针变量将会占用多出一倍的寄存器空间。由于Tesla C2050的设备存储器容量不超过4 GB（在增加ECC的情况下可用的设备存储器空间只有2.625 GB），因此完全没有必要在设备代码中使用64 b的指针。因此，在本文实现的Fermi体系结构上的GEMM代码中，CPU代码和设备代码总是被分别编译的。

1.4 设备存储器访问

在GT200体系结构上，对设备存储器的访问是按照半个warp的单位来进行处理的，而在Fermi体系结构上，是按照一个warp的单位来进行处理的。因此，程序员需要调整内核调用时的维度设置。对于具有两个维度的线程块，其x维度的大小应该是warp大小的整倍数，而非半个warp大小的整倍数，从而使得每一个warp在访问设备存储器时可以得到较高的性能。

1.5 Bank冲突

在GT200体系结构上，共享存储器具有16个bank，而且对其的访问是按照半个warp的单位来进行处理的，而在Fermi体系结构上，共享存储器具有32个bank，而且对其的访问是按照整个warp的单位来进行处理的。每一个bank的大小为32位。

1.6 多内核并行执行

Fermi体系结构支持多个内核代码的并行执行，使得不同应用上下文的内核代码可以同时在一个GPU上运行；这样，多个小的内核代码可以共同利用一个GPU上的计算资源。这也是Fermi体系结构的新特性，但在本文自动优化的GEMM代码中并未使用到。

2 自动优化的GEMM代码

2.1 GEMM代码模板

文献[4]描述了在GTX280实现的达到393 Gflops性能的SGEMM内核代码，在此，依然使用该代码作为实现自动优化的代码模板。在该实现中，一个Csub被一个线程块进行计算，根据线程块中线程数目的多少，一个线程可以计算Csub的半列或多列元素。例如，对于m= 16和n=64，而线程块具有16×4个线程，则每一个线程将计算一整列Csub的元素；如果线程块具有16×8个线程，则每一个线程将计算半列Csub的元素。与文献[4]中的实现类似，每一个线程在使用线程ID计算出访问矩阵A、B和C的指针位置之后，进入一个循环。在每一轮循环中，一个线程块从设备存储器中读入一个Asub的数据到共享存储器中，之后，又一个内层的循环被执行：在每一轮内层循环中，一个线程从矩阵B中读入一个或多个元素，把这些元素与共享存储器中相应的数据做计算，将结果累加到Csub对应的寄存器中。最后，每个线程将其计算的Csub的数据写回到设备存储器中。

尽管与文献[4]中的GEMM实现具有类似的代码结构，但是为了满足缓存的友好性，对该代码模板进行了一个重要的修正。在此计算中，矩阵B总是位于设备存储器中，一列线程块总是需要读取一列Bsub的数据，如图2（a）所示。如果线程块能够按照列优先的顺序调度，则可以达到更好的缓存命中效果。因此，在此GEMM模板中，blockIdx.x和blockIdx.y变量的顺序被对调，从而达到如图2（b）所示的设备存储器访问效果。这样，在进行这个转置之后，线程块被以较优的缓存命中效果的方式调度。在后面给出的性能测试中，将对比这个对调进行或不进行的性能结果。

2.2 对代码模板进行自动优化

设计的自动调节程序根据代码模板生成代码并测试其性能结果。在代码模板中，5个参数（m，k，n，tx和ty）决定代码模板的执行行为；此外，引入2个额外的参数，一个用来决定线程块的维度是否对调，来测试缓存的效果，另一个确定L1缓存和共享存储器的配置比例。因此，整个代码模板的行为是由7个参数确定的。

图3显示经过优化得到的矩阵大小为2 048时的DGEMM代码的实例。注意：在该代码中，blockIdx.x和blockIdx.y变量通过C语言的宏进行了对调，而相应的调用代码也做了相应的改变。在此代码中，m=8，k= 64，n=1 024，tx=64，ty=8，因此一个线程负责计算Csub中的两列数据。

3 性能测试

4 结 语

为了在Fermi体系结构上书写高效的代码，程序员需要很好地了解Fermi体系结构的新特性，以及它们是如何影响程序的性能的。对于Fermi体系结构，程序员尤其要关心缓存对性能的影响。对于Fermi这种复杂的行为和性能难以预测的硬件体系结构，自动优化技术不失为得到高性能代码的一种实用技术。

参考文献

[1] NVIDIA Corp. Whitepaper： NVIDIA's next generation CUDA compute architecture [R/OL]. [2012?05?18]. http：//www. insidehpc.com.

[2] NVIDIA Corp. Tuning CUDA applications for Fermi [R/OL]. [2011?05?03]. http：// www. people.maths.ox.ac.uk/gilesm/cuda/doc/Fermi_Tuning_Guide.

[3] NVIDIA Corp. CUDA compute unified device architecture， programming guide， Version 3.0 [R/OL]. [2010?05?03]. http：// www. mohamedfahmed.wordpress.com

[4] CUI Xiang， CHEN Yi?feng， MEI Hong， et al. Auto?tuning GEMM for GPGPU with Cache [C]// Proceedings of 2010 IEEE 16th International Conference on Parallel and Distributed Systems （ICPADS）. Shanghai， China： IEEE， 2010： 237?242.

[5] VOLKOV V， DEMMEL J W. Benchmarking GPUs to tune dense linear algebra [C]// Proceedings of 2008. International Conference for High Performance Computing， Networking， Storage and Analysis Austin， TX： [s.n.]， 2008： 1?11.

[6] 李晓雯，崔翔.GPU矩阵乘法和FFT算法的性能优化[J].现代电子技术，2013，36（4）：80?84.

1.3 32/64位设备代码

在Fermi体系结构上，如果代码按照64?bit的模式编译，则CUDA编译器会将CPU代码和设备代码都编译成为64?bit的目标代码。在这种情况下，设备代码中的指针变量将会占用多出一倍的寄存器空间。由于Tesla C2050的设备存储器容量不超过4 GB（在增加ECC的情况下可用的设备存储器空间只有2.625 GB），因此完全没有必要在设备代码中使用64 b的指针。因此，在本文实现的Fermi体系结构上的GEMM代码中，CPU代码和设备代码总是被分别编译的。

1.4 设备存储器访问

在GT200体系结构上，对设备存储器的访问是按照半个warp的单位来进行处理的，而在Fermi体系结构上，是按照一个warp的单位来进行处理的。因此，程序员需要调整内核调用时的维度设置。对于具有两个维度的线程块，其x维度的大小应该是warp大小的整倍数，而非半个warp大小的整倍数，从而使得每一个warp在访问设备存储器时可以得到较高的性能。

1.5 Bank冲突

在GT200体系结构上，共享存储器具有16个bank，而且对其的访问是按照半个warp的单位来进行处理的，而在Fermi体系结构上，共享存储器具有32个bank，而且对其的访问是按照整个warp的单位来进行处理的。每一个bank的大小为32位。

1.6 多内核并行执行

Fermi体系结构支持多个内核代码的并行执行，使得不同应用上下文的内核代码可以同时在一个GPU上运行；这样，多个小的内核代码可以共同利用一个GPU上的计算资源。这也是Fermi体系结构的新特性，但在本文自动优化的GEMM代码中并未使用到。

2 自动优化的GEMM代码

2.1 GEMM代码模板

文献[4]描述了在GTX280实现的达到393 Gflops性能的SGEMM内核代码，在此，依然使用该代码作为实现自动优化的代码模板。在该实现中，一个Csub被一个线程块进行计算，根据线程块中线程数目的多少，一个线程可以计算Csub的半列或多列元素。例如，对于m= 16和n=64，而线程块具有16×4个线程，则每一个线程将计算一整列Csub的元素；如果线程块具有16×8个线程，则每一个线程将计算半列Csub的元素。与文献[4]中的实现类似，每一个线程在使用线程ID计算出访问矩阵A、B和C的指针位置之后，进入一个循环。在每一轮循环中，一个线程块从设备存储器中读入一个Asub的数据到共享存储器中，之后，又一个内层的循环被执行：在每一轮内层循环中，一个线程从矩阵B中读入一个或多个元素，把这些元素与共享存储器中相应的数据做计算，将结果累加到Csub对应的寄存器中。最后，每个线程将其计算的Csub的数据写回到设备存储器中。

尽管与文献[4]中的GEMM实现具有类似的代码结构，但是为了满足缓存的友好性，对该代码模板进行了一个重要的修正。在此计算中，矩阵B总是位于设备存储器中，一列线程块总是需要读取一列Bsub的数据，如图2（a）所示。如果线程块能够按照列优先的顺序调度，则可以达到更好的缓存命中效果。因此，在此GEMM模板中，blockIdx.x和blockIdx.y变量的顺序被对调，从而达到如图2（b）所示的设备存储器访问效果。这样，在进行这个转置之后，线程块被以较优的缓存命中效果的方式调度。在后面给出的性能测试中，将对比这个对调进行或不进行的性能结果。

2.2 对代码模板进行自动优化

设计的自动调节程序根据代码模板生成代码并测试其性能结果。在代码模板中，5个参数（m，k，n，tx和ty）决定代码模板的执行行为；此外，引入2个额外的参数，一个用来决定线程块的维度是否对调，来测试缓存的效果，另一个确定L1缓存和共享存储器的配置比例。因此，整个代码模板的行为是由7个参数确定的。

图3显示经过优化得到的矩阵大小为2 048时的DGEMM代码的实例。注意：在该代码中，blockIdx.x和blockIdx.y变量通过C语言的宏进行了对调，而相应的调用代码也做了相应的改变。在此代码中，m=8，k= 64，n=1 024，tx=64，ty=8，因此一个线程负责计算Csub中的两列数据。

3 性能测试

4 结 语

为了在Fermi体系结构上书写高效的代码，程序员需要很好地了解Fermi体系结构的新特性，以及它们是如何影响程序的性能的。对于Fermi体系结构，程序员尤其要关心缓存对性能的影响。对于Fermi这种复杂的行为和性能难以预测的硬件体系结构，自动优化技术不失为得到高性能代码的一种实用技术。

参考文献

[1] NVIDIA Corp. Whitepaper： NVIDIA's next generation CUDA compute architecture [R/OL]. [2012?05?18]. http：//www. insidehpc.com.

[2] NVIDIA Corp. Tuning CUDA applications for Fermi [R/OL]. [2011?05?03]. http：// www. people.maths.ox.ac.uk/gilesm/cuda/doc/Fermi_Tuning_Guide.

[3] NVIDIA Corp. CUDA compute unified device architecture， programming guide， Version 3.0 [R/OL]. [2010?05?03]. http：// www. mohamedfahmed.wordpress.com

[4] CUI Xiang， CHEN Yi?feng， MEI Hong， et al. Auto?tuning GEMM for GPGPU with Cache [C]// Proceedings of 2010 IEEE 16th International Conference on Parallel and Distributed Systems （ICPADS）. Shanghai， China： IEEE， 2010： 237?242.

[5] VOLKOV V， DEMMEL J W. Benchmarking GPUs to tune dense linear algebra [C]// Proceedings of 2008. International Conference for High Performance Computing， Networking， Storage and Analysis Austin， TX： [s.n.]， 2008： 1?11.

[6] 李晓雯，崔翔.GPU矩阵乘法和FFT算法的性能优化[J].现代电子技术，2013，36（4）：80?84.