GPU矩陣乘法和FFT算法的性能優(yōu)化

摘 要： 當前GPU的體系結構為高性能計算提供了良好的可編程性。為了得到眾核GPU高性能程序設計的一般方法，探索GPU程序性能優(yōu)化技術，對在GPU上進行高性能程序設計的經(jīng)驗進行了總結。通過基準測試，得到GPU性能指標，對GPU程序設計進行指導。使用CUDA對單精度矩陣乘法和FFT進行性能優(yōu)化，前一個算法是計算密集型任務，后一個算法是帶寬密集型任務。在NVIDIA GeForce GTX280 GPU上，矩陣乘法算法達到393 Gflop/s的峰值速度，比CUBLAS 2.0數(shù)學庫提高了5%；對于一些維度的FFT計算也取得了較好的性能。

關鍵詞： GPU程序設計； 矩陣乘法； 快速傅里葉變換； 性能優(yōu)化技術

中圖分類號： TN911?34； TP312 文獻標識碼： A 文章編號： 1004?373X（2013）04?0080?05

0 引 言

當前GPU的體系結構為高性能計算提供了良好的可編程性[1?6]。一個GPU包含多個多處理器（MP），一個多處理器又包含若干個流處理器（SP）。在GPU上具有全局設備存儲器可以被所有多處理器訪問；每一個多處理器又具有一定量的共享存儲器（Shared Memory）供流處理器共享使用和一定量的寄存器供流處理器獨占使用。CUDA編程模型為這樣的硬件提供了易用的編程方式[7]。它的運行模型包括若干不同的并行層次：線程塊網(wǎng)格（Thread Block Gird）、線程塊（Thread Block）、Warp塊以及線程，它們被調度到GPU上執(zhí)行。一個線程塊只會被調度到一個多處理器上執(zhí)行，因此，一個線程塊中的線程可以共享訪問一部分共享存儲器。一個線程塊網(wǎng)格包括很多線程塊。因為線程塊的數(shù)量可以大大超過多處理器的數(shù)量，一個線程塊網(wǎng)格中的線程塊被動態(tài)的調度到多處理器上運行。線程塊中的線程被組織成Warp塊，一個Warp塊以SIMD的方式執(zhí)行。線程塊中的線程以排他的方式被調度到多處理器上執(zhí)行。

CUDA提供了對眾核GPU進行程序設計的底層方式，但是要想得到頂級的性能對程序員來講確實是一個挑戰(zhàn)。程序員必須理解GPU硬件平臺的各種制約因素和限制來設計算法，并在影響性能的各種因素之間進行權衡。例如，早期硬件生產(chǎn)商提供的GPU上的數(shù)學庫：

（1）基本線性代數(shù)運算的CUBLAS庫；

（2）FFT計算的CUFFT庫。

其中的一些基本運算函數(shù)（版本為CUBLAS 1.0和CUFFT 1.1）就具有比較差的性能。在本文中，對兩個基本的運算函數(shù)：計算密集型的矩陣乘法和帶寬/通信密集型的1維FFT進行性能優(yōu)化，并以此說明GPU程序優(yōu)化技術。

本文優(yōu)化的代碼的性能遠遠優(yōu)于早期硬件生產(chǎn)商提供的GPU上的數(shù)學庫。這說明，對于一個程序員而言，使用CUDA在GPU上得到頂級性能的代碼是十分困難的。實際上，程序員不僅要精確了解GPU的硬件體系結構，還要具有豐富的實踐經(jīng)驗，來幫助他們達到頂級的存儲器帶寬、峰值速度等性能指標。在本文中，通過一些基本的測試，得到了相應的編程原則。本文的代碼性能為單精度浮點數(shù)在NVIDIA GeForce GTX280上得到的性能。

1 基本優(yōu)化原則

一系列的GTX200 GPU程序設計優(yōu)化原則可以在文獻[2?4]中找到。在本節(jié)中，主要討論根據(jù)所進行的基礎實驗，得到的一些優(yōu)化原則。

1.1 以MIMD方式控制線程Warp

在每一個時鐘周期，一個多處理器以SIMD的方式運行：即一個多處理器中的每一個微處理器執(zhí)行相同的指令運算不同的數(shù)據(jù)[8]。基本的SIMD的單位是Warp，每一個Warp具有相同的線程個數(shù)。在GTX200 GPU上，Warp的大小是32。Warp塊是按照線程的自然順序形成的：前32個線程是第1個Warp，而后第2個等。活動的Warp以時間片的方式占用多處理器的計算資源。

在大多數(shù)的CUDA程序中，每一個線程塊中的線程具有相同的代碼結構，只是根據(jù)線程ID區(qū)分處理不同的數(shù)據(jù)元素。盡管沒有提倡，但實際上不同的Warp塊可以根據(jù)不同的Warp ID執(zhí)行不同的代碼。這通過在條件語句中測試Warp ID來實現(xiàn)。只要同一個Warp塊中的線程執(zhí)行同樣的條件分支，在性能上就不會有明顯的損失。這使得代碼可以運行在MIMD的模式下并執(zhí)行不同的指令。

圖1演示了這種MIMD機制的實現(xiàn)方式。在這種情況下，使用這種機制可以減少一個線程塊的共享存儲器的使用量，從而使得更多的線程塊可以被加載到一個多處理器上。在FFT計算中，各個線程需要進行全局的矩陣轉置：4個Warp塊中的各個線程將自己的數(shù)據(jù)放入到共享存儲器中并同步，之后再從共享存儲器中將自己需要的數(shù)據(jù)取回，如圖1（a）所示。如果按照圖1（b）所示的方式進行，共享存儲器的用量可以被減半：前2個Warp塊中的線程先將自己的數(shù)據(jù)放入到共享存儲器中并同步，各個線程將自己需要的數(shù)據(jù)取回，而后后2個Warp塊中的線程再將自己的數(shù)據(jù)放入到共享存儲器中并同步，各個線程將自己需要的數(shù)據(jù)取回。根據(jù)實驗，具有Warp間條件指令轉移的代碼會帶來性能10%～30%的損失。但是，由于一個多處理器可以具有更多的活動的線程塊，在本文FFT的代碼實現(xiàn)中，損失的性能可以被由于具有更多的活動線程塊掩蓋指令流水線延遲得到的性能所掩蓋。

1.2 峰值指令吞吐量

通過測試代碼可以測試一個GPU硬件所能達到的峰值指令吞吐量，通過比較該值和代碼達到的性能之間的差異，可以推測該代碼進一步優(yōu)化的空間。GeForce GTX280具有240個微處理器，運行在1.295 GHz的頻率下，單精度浮點數(shù)的理論性能為933 Gflop/s。通過在一個展開的循環(huán)中運行寄存器間的MAD指令a=b*a+b和b=a*b+a，可以達到GPU的峰值指令吞吐量。在完全流水線的情況下，可以達到622 Gflop/s的性能。

1.3 流水線延遲

在文獻[2?4]中指出，每一個多處理器必須要有192個活動的線程，才能夠掩蓋流水線的延遲。本文的實驗表明，必須至少要有256個活動的線程才能夠掩蓋流水線的延遲。只有在每一個多處理器有超過256個活動線程，而不論這256個活動線程分成多少個線程塊，峰值指令吞吐量620 Gflop/s才可能被達到。

1.4 增加活動線程塊以掩蓋同步指令

表1示意了不同同步和計算指令比例以及在一個多處理器上具有1個或2個活動線程塊的情況下，代碼可以達到的指令吞吐量。如果一個多處理器上只有一個活動的線程塊，代碼中的同步指令會導致一個多處理器成為空閑的狀態(tài)。如果一個多處理器上有多個活動的線程塊，在一個線程塊執(zhí)行同步指令的時候，另一個活動的線程塊可以被調度器自動調度執(zhí)行。這種使用多個活動線程塊掩蓋線程塊同步指令的方法也被運用于FFT的代碼優(yōu)化中。

表1 掩蓋同步指令延遲

1.5 設備存儲器帶寬

設備存儲器不具有緩存，所以必須使用凝聚的訪問方式以得到最大的設備存儲器帶寬[8]。GPU可以使用一個指令將32 b，64 b或128 b的字由設備存儲器讀入到寄存器中，而半個warp塊中的線程必須以凝聚的訪問方式讀寫設備存儲器。

通過測試設備存儲器峰值帶寬，凝聚的64 b的訪問性能略強于32 b字的訪問，而在每個線程塊具有128個線程的情況下，32 b字和64 b字訪問的峰值帶寬分別是121.6 GB/s和124.8 GB/s（注意增加_syncthreads（）指令可以略微增加帶寬性能）。在矩陣乘法的實現(xiàn)中采用每個線程塊具有128個線程的設計，這樣基本上可以保證設備存儲器訪問的最高性能。

2 SGEMM實現(xiàn)和性能測試

上述程序優(yōu)化原則可以被用于單精度矩陣乘法（SGEMM）實現(xiàn)。SGEMM屬于計算密集型的任務，但在設計不當?shù)那闆r下，會受到設備存儲器帶寬的約束。

計算矩陣A和矩陣B的乘積得到矩陣C的計算可以按照下述分塊的方式進行：每一個線程塊負責計算矩陣C的一個子塊Csub，而該線程塊中的每一個線程負責計算子塊Csub中的若干個元素。Csub是兩個長方形矩陣的乘積：矩陣A的和Csub具有相同行數(shù)的長方形矩陣乘以矩陣B的和Csub具有相同列數(shù)的長方形矩陣，如圖2所示。

由于GPU硬件存儲的限制，這兩個長方形的矩陣也被分割為相應的子塊，稱為Asub和Bsub，Csub為對應Asub和Bsub子塊矩陣乘積的和。Csub的大小決定了每輪運算必須被加載到GPU芯片的數(shù)據(jù)量大小。如果Csub過小，則SGEMM會變成一個帶寬受限的問題。

假設Asub，Bsub和Csub的大小分別為m×k，k×n和m×n，而矩陣A，B和C分別被分為M×K，K×N和M×N個這樣的子塊。計算一個Csub需要從設備存儲器中取K個Asub和Bsub。由于在矩陣C中有M×N個Csub，則整個計算需要從設備存儲器中取M×m×N×n×K×k×（[1m]+[1n]）個浮點數(shù)。

該算法以行優(yōu)先的方式設計。線程塊的線程數(shù)目為128，被組織成4×32的格式。Csub的大小為16×256。每一個線程塊負責計算一個Csub，而一個線程負責計算Csub中的2列。圖3顯示了每一個線程的流程圖。

首先，線程塊申請16×32的共享存儲器空間用于存儲矩陣A的子塊，每一個線程申請寄存器用于存儲Csub中的16個float2即32個float。在使用線程ID計算出訪問矩陣A，B和C的指針位置之后，進入一個循環(huán)。在每一輪循環(huán)中，一個線程塊從設備存儲器中讀入16×32的Asub的數(shù)據(jù)到共享存儲器中，之后，又一個內層的循環(huán)被執(zhí)行：在每一輪內層循環(huán)中，一個線程從矩陣B中讀入一個float2，把這兩個float和共享存儲器中相應的數(shù)據(jù)做計算，將結果累加到Csub對應的寄存器中。最后，每個線程將其計算的2列Csub的數(shù)據(jù)寫回到設備存儲器中。

圖4對本文的代碼和CUBLAS代碼的性能進行了比較。當矩陣大小為8 192×8 192（在GTX285 GPU上可以進行計算的最大尺寸）時，本文的代碼達到了393 Gflop/s的性能，比文獻[8]中實現(xiàn)的性能大約高5%。而使用共享存儲器進行MAD計算的理論指令峰值吞吐量為410 Gflop/s。本文的實現(xiàn)和文獻[8]中的SGEMM實現(xiàn)具有類似的代碼結構，下列因素導致本文的代碼取得了更優(yōu)的性能：

（1）在本文的代碼中，Csub的尺寸為16×256，而非16×64，從而減少了從設備存儲器訪問數(shù)據(jù)的數(shù)量；

（2）在本文的代碼中，一個線程塊包含128個線程，而每一個線程以凝聚方式訪問一個float2即64 b字，這樣可以取得較高的設備存儲器帶寬；

（3）在本文的代碼中，一個Asub的大小為16×32，而非16×16，從而使得整個算法的同步指令次數(shù)被大大減小了。

3 FFT實現(xiàn)和性能測試

另一個代碼優(yōu)化的示例為批量的1維復數(shù)FFT，根據(jù)向量的大小不同，1維復數(shù)FFT的實現(xiàn)可以被分為三種情況：當向量的長度為8或16時，不需要進行線程之間的數(shù)據(jù)交換，一個向量的計算可以在寄存器中被一個線程完成；當向量的長度為64，256，512，1 024，2 048或4 096時，線程之間需要進行數(shù)據(jù)交換，且數(shù)據(jù)交換可以在共享存儲器中完成；當向量的長度超過8 192時，必須通過設備存儲器進行矩陣數(shù)據(jù)的轉置，這可以通過調用兩次內核代碼完成。因此，當向量的長度為8或16時，此計算為存儲帶寬受限型；當向量的長度為2 048或4 096時，此計算為通信受限型。在此，取向量的長度為16，2 048或4 096兩種情況來說明程序優(yōu)化的方法。

3.1 向量長度為16

當向量長度為16時，這是一個存儲帶寬受限型的問題。為了得到良好的性能，必須要使得對設備存儲器的訪問達到一個良好的帶寬。在進行計算之前，原始數(shù)據(jù)被首先進行整理，使得線程可以以凝聚的方式，讀取一個向量中需要計算的16個元素。盡管這種數(shù)據(jù)整理可以通過在共享存儲器中的數(shù)據(jù)轉置來代替，但在本文不涉及到此問題。所以在這個實現(xiàn)中，一個線程需要做的就是從設備存儲器中讀取16個元素，進行1維FFT計算，再把數(shù)據(jù)寫回設備存儲器。由于每一個元素為復數(shù)，所以每一個線程自然的就在讀取64 b字的數(shù)據(jù)。此外，根據(jù)前面提到的優(yōu)化原則，在每一個線程塊中包含256個線程，并且在計算和讀寫設備存儲器之間加入同步指令，以得到最優(yōu)的設備存儲器訪問帶寬。

3.2 向量長度為2 048或4 096

當向量長度為2 048時，每一個線程塊包含128個線程，每一個多處理器上可以具有2個活動線程塊。每一個線程負責計算2 048個復數(shù)元素，這通過在寄存器中計算長度為16，16和8的向量完成，在3次計算之間，通過共享存儲器交換數(shù)據(jù)。為了在共享存儲器中交換2 048個單精度復數(shù)數(shù)據(jù)，至少需要8 KB的空間用來交換復數(shù)的實部或虛部，使得一個多處理器上只能存在1個活動線程塊。前面的優(yōu)化原則提到，一個多處理器上至少需要2個活動線程塊以掩蓋同步指令的開銷。本文的解決方法是，通過以下操作將共享存儲器的使用量減少一半：首先前一半的線程先將自己的數(shù)據(jù)放入到共享存儲器中并同步，各個線程將自己需要的數(shù)據(jù)取回，而后后一半的線程再將自己的數(shù)據(jù)放入到共享存儲器中并同步，各個線程將自己需要的數(shù)據(jù)取回。這樣，一個多處理器上可以同時具有兩個活動的線程塊。在此實現(xiàn)中，通過控制Warp ID來達到此種操作效果。當向量長度為4 096時，實現(xiàn)的思路與向量長度為2 048是類似的，只不過此時一個多處理器上只能具有一個活動的線程塊。

3.3 性能測試

對64 MB的數(shù)據(jù)進行批量1維FFT的性能測試：當向量大小為1 024時，對應的信號數(shù)目為8 129；當向量大小為2 048時，對應的信號數(shù)目為4 096，等。圖5比較了本文的代碼和CUFFT 1.1，Volkov和Kazian的代碼[5?6]的性能。注意，在本文的實現(xiàn)中，當向量大小為2 048和4 096時，該計算通過一次內核調用完成，從而使得數(shù)據(jù)不需要進行額外的整理，而其他一些實現(xiàn)則要求數(shù)據(jù)在GPU計算完成之后使用CPU在主存中對數(shù)據(jù)進行整理。

4 結 語

書寫高性能的GPU代碼是十分困難的，本文介紹的工作并不是為了得到特定算法的高性能代碼，而是為了總結幫助程序員得到高性能GPU代碼的程序設計原則。程序員在進行算法設計時，往往需要進行自頂向下和自下向上的反復優(yōu)化。本文中提到的優(yōu)化原則可以作為對CUDA程序設計指南和其他一些工作的補充。在硬件體系結構發(fā)生變化的情況下，原先非最優(yōu)的代碼可能又會取得頂級的性能。取得特定算法的高性能代碼本身并不是本文所看重的，而在進行代碼優(yōu)化過程中取得的經(jīng)驗，有助于對GPU和GPU集群程序設計問題進行進一步的研究。

參考文獻

[1] 吳恩華，柳有權.基于圖形處理器（GPU）的通用計算[J].計算機輔助設計與圖形學學報，2004，16（5）：601?612.

[2] 吳恩華.圖形處理器用于通用計算的技術、現(xiàn)狀及其挑戰(zhàn)[J].軟件學報，2004，15（10）：101?112.

[3] 蘇超軾，趙明昌，張向文.GPU加速的八叉樹體繪制加速算法[J].計算機應用，2008（10）：51?62.

[4] 儲璟駿，楊新，高艷.使用GPU編程的光線投射體繪制算法[J].計算機輔助設計與圖形學學報，2007，19（5）：257?262.

[5] 朱亞平，楊慧珠，董淵，等.OpenGL技術在地震數(shù)據(jù)可視化中的應用[J].石油地球物理勘探，2000（8）：35?42.

[6] 徐品，藍善禎，劉蘭蘭.利用GPU進行通用數(shù)值計算的研究[J].中國傳媒大學學報：自然科學版，2009（4）：101?112.

[7] NVIDIA Corporation. NV1DIA CUDA programming guide Version 0.1 [EB/OL]. [2009?04?28]. http：// developer.nvidia.con/cuda.

[8] COPPERSMITH D， WINOGRAD S. Matrix multipli cation via arithmetic progressions [C]// Proceedings of the nineteenth annual ACM symposium on Theory of computing. New York， NY， USA： ACM， 1987： 1?6.

[9] RYOO S. Optimization principles and application performance evaluation of a multithreaded GPU using CUDA [C]// Proceedings of 13th ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming. [S.l.]： ACM， 2008： 73?82.

[10] VOLKOV V， DEMMEL J W. Benchmarking GPUs to tune dense linear algebra SC08 [C]// Proceedings of the 2008 ACM/IEEE Conference on Supercomputing. [S.l.]： IEEE Press， 2008： 1?11.

[11] GOTO K， VAN DE GEIJN R A. Anatomy of high performance matrix multiplication [J]. ACM Transactions on Mathematical Software， 2008， 3（43）： 12?25.

[12] HWU W M， RODRIGUES C， RYOO S， et al. Compute unified device architecture application suitability [J]. Computing in Science and Engineering， 2009， 11（3）： 16?26.