摘  要： 首先介紹了CUDA架構(gòu)特點，在GPU上基于CUDA使用兩種方法實現(xiàn)了矩陣乘法，并根據(jù)CUDA特有的軟硬件架構(gòu)對矩陣乘法進行了優(yōu)化。然后計算GPU峰值比并進行了分析。實驗結(jié)果表明，基于CUDA的矩陣乘法相對于CPU矩陣乘法獲得了很高的加速比，最高加速比達到1 079.64。GPU浮點運算能力得到有效利用，峰值比最高達到30.85%。
關鍵詞： CUDA；矩陣乘法；加速比；峰值比

    隨著多核CPU和眾核GPU的快速發(fā)展，計算行業(yè)正在從只使用CPU的“中央處理”向CPU與GPU并用的“協(xié)同處理”發(fā)展，并行系統(tǒng)已成為主流處理器芯片。傳統(tǒng)的GPU架構(gòu)受其硬件架構(gòu)的影響不能有效利用其資源進行通用計算，NVIDIA(英偉達)公司推出的統(tǒng)一計算設備架構(gòu)CUDA(Compute Unified Device Architecturem)，使得GPU具備更強的可編程性，更精確和更高的性能，應用領域也更加廣泛。
    矩陣乘法是一種大計算量的算法，也是很耗時的運算。CPU提高單個核心性能的主要手段比如提高處理器工作頻率及增加指令級并行都遇到了瓶頸，當遇到運算量大的計算，CPU進行大矩陣的乘法就變得相當耗時,運算效率很低下。因此，GPU憑借其超強計算能力應運而生，讓個人PC擁有了大型計算機才具備的運算能力。本文運用GPU的超強計算能力在CUDA架構(gòu)上實現(xiàn)了大矩陣乘法。
1 CUDA架構(gòu)
    NVIDIA及時推出CUDA這一編程模型，在應用程序中充分結(jié)合利用CPU和GPU各自的優(yōu)點，特別是GPU強大的浮點計算能力。CPU主要專注于數(shù)據(jù)高速緩存(cache)和流處理(flow control)，而GPU更多地專注于計算密集型和高度并行的計算。盡管GPU的運行頻率低于CPU，但GPU憑著更多的執(zhí)行單元數(shù)量使其在浮點計算能力上獲得較大優(yōu)勢[1]。當前的NVIDIA GPU中包含完整前端的流多處理器(SM)，每個SM可以看成一個包含8個1D流處理器(SP)的SIMD處理器。主流GPU的性能可以達到同期主流CPU性能的10倍左右。圖1所示為GPU與CPU峰值浮點計算能力的比較。

    CUDA的編程模型是CPU與GPU協(xié)同工作，CPU作為主機(Host)主要負責邏輯性強的事務處理及串行計算，GPU作為協(xié)處理器或者設備(Device)負責密集型的大規(guī)模數(shù)據(jù)并行計算。一個完整的CUDA程序=CPU串行處理+GPU Kernel函數(shù)并行處理。
    一個CUDA架構(gòu)下的程序分為兩個部分，即上述的Host端和Device端。通常情況下程序的執(zhí)行順序如下：Host端程序先在CPU上準備數(shù)據(jù)，然后把數(shù)據(jù)復制到顯存中，再由GPU執(zhí)行Device端程序來處理這些數(shù)據(jù)，最后Host端程序再把結(jié)束運算后的數(shù)據(jù)從顯存中取回。
    圖2為CUDA編程模型，從中可以看出，Thread是GPU執(zhí)行運算時的最小單位。也就是說，一個Kernel以線程網(wǎng)格Grid的形式組織，每個Grid由若干個線程塊Block組成，而每個線程塊又由若干個線程Thread組成。一個Kernel函數(shù)中會存在兩個層次的并行，Grid中Block之間的并行和Block中Thread之間的并行，這樣的設計克服了傳統(tǒng)GPGPU不能實現(xiàn)線程間通信的缺點[2]。

    同一個Block下的Thread共用相同的共享存儲器，通過共享存儲器交換數(shù)據(jù)，并通過柵欄同步保證線程間能夠正確地共享數(shù)據(jù)。因此，一個Block下的Thread雖然是并行的，但在同一時刻執(zhí)行的指令并不一定都相同，實現(xiàn)了不同Thread間的協(xié)同合作。這一特性可以顯著提高程序的執(zhí)行效率，并大大拓展GPU的適用范圍。
2 基于CUDA架構(gòu)矩陣乘法的實現(xiàn)
2.1 一維帶狀劃分
    給定一個M×K的矩陣A和一個K×N的矩陣B，將矩陣B乘以矩陣A的結(jié)果存儲在一個M×N的矩陣C中。此種矩陣乘法使用了一維帶狀劃分，每個線程將負責讀取矩陣A中的一行和B中的一列，矩陣進行乘法運算并將計算結(jié)果存儲在全局存儲器。
    全局存儲器會對矩陣A進行N次讀取，對矩陣B進行M次讀取。假設數(shù)組在每個維度上的尺寸都是BLOCK_SIZE的整數(shù)倍。若矩陣大小為32×32，則可表示為(2×16)×(2×16)。下面的內(nèi)核定義中，結(jié)果矩陣C中的每個元素由一個線程負責，for()循環(huán)完成A中第X行與B中第X列對應元素的乘加運算，并將結(jié)果累加到Cvalue。
      For( int e=0;e < A.width;++e)
      Cvalue+=A.elements[row*A.width+e] *
B.elements[e*B.width+col];
      C.elements[row*width+col]=Cvalue;
      在矩陣相乘實現(xiàn)中，這個內(nèi)核運算的速度不盡人意，主要瓶頸在于對內(nèi)存的重復讀取，計算量是2×M×N×K flop，而全局內(nèi)存的訪問量為2×M×N×K B[3]。若矩陣維數(shù)為1 024×1 024，則此次矩陣相乘的計算量就有2 G flop，當矩陣維數(shù)更大時，這個運算量就相當大，在內(nèi)存的讀取上會浪費大量的時間。
2.2 二維棋盤劃分
    因為矩陣A的行和矩陣B的列多次被讀取，為了避免重復加載，選擇把矩陣進行分塊運算，使用shared memory來實現(xiàn)矩陣乘法。運用shared memory的好處在于其延遲小于global memory，并且還能使線程間進行通信。矩陣A只被讀了N/BLOCK_SIZE次，矩陣B僅被讀了M/BLOCK_SIZE次，節(jié)省了大量的global memory帶寬。
    首先把劃分的小矩陣塊加載到share memory，則小矩陣本身的乘法就不用去存取外部的任何內(nèi)存了，因此在二維棋盤劃分中，矩陣乘法的計算量仍然是2×M×N×K flop，b是矩陣B劃分的小矩陣塊的大小，則全局內(nèi)存訪問量是2×M×N×K/b B。
    棋盤劃分運算可以表示為：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)+…+A(0～15,(16×(n-1)-1)～(16×(n-1))×B((16×(n-1)-1)～(16×(n-1)),0～15)。
        for (int j=0;j<wA;j+=BLOCK_SIZE)
    {    //聲明用于存儲A，B子塊的share memory數(shù)組
        __shared__ float As[BLOCK_SIZE][BLOCK_SIZE];
        __shared__ float Bs[BLOCK_SIZE][BLOCK_SIZE];
        }...
        //兩個子塊的乘加，每個線程負責C中一個元素
值的計算
        for (int k = 0; k < BLOCK_SIZE; ++k)
        {
           float t;
           C sub + = As [ ty ] [ k ] * Bs [ k ] [ tx ];
           Cs [ ty ] [ tx ] = C sub;
        }
    __syncthreads();
    ....
    C[(by*BLOCK_SIZE+ty)*wA+bx*BLOCK_SIZE+tx] = Csub;
    dim3 myblock(BLOCK_SIZE,BLOCK_SIZE,1);
    dim3mygrid(((wB+BLOCK_SIZE-1)/BLOCK_SIZE),
(wB+BLOCK_SIZE-1)/BLOCK_SIZE,1);
    根據(jù)NVIDIA CUDA Programming Guide，一個Block 里至少要有64個Thread，最多有512個Thread。官方建議256個Thread是最合適的，因為此時有足夠多的active warp有效地隱藏延遲，使得SM能夠盡量滿負荷工作[4]。為便于理解，假設矩陣為n×n ,此時BLOCK_SIZE設置為16，使用dim3來設計，每個Block包含16×16個Thread，一個Grid共有(n/16)×(n/16)個Block。
    BLOCK_SIZE是不是越大越好呢？這樣一個SM里的Thread 就更多，雖然Thread越多越能隱藏 latency，但G80/G92架構(gòu)每個SM上shared memory僅有16 KB，這會讓每個Thread 能使用的資源更少，效率反而會下降。
2.3 根據(jù)CUDA架構(gòu)對矩陣乘法進行優(yōu)化

    因為棋盤劃分中涉及到的是二維數(shù)組，cudaMalloc2D()能確保分配二維數(shù)組并且能分配適當?shù)奶畛湟詽M足對齊要求，還能確保在訪問行地址或者二維數(shù)組與其他設備內(nèi)存之間的數(shù)據(jù)復制能達到最佳性能。
    二維棋盤劃分方法僅限于數(shù)組大小必須是BLOCK_SIZE的整數(shù)倍，若矩陣維數(shù)并不是16的整數(shù)倍，則會造成運算效率的下降，此時可以利用CUDA架構(gòu)特點和CUDA提供的cudaMallocPitch()函數(shù)來解決此問題。cudaMallocPitch()可以自動地以最佳倍數(shù)來分配內(nèi)存。
    呼叫Kernel部分需要修改成：
    matrixMul<<<mygrid,myblock>>>(d_A,d_B,d_C,wA,wB,
d_pitchA/sizeof(float),d_pitchB/siz
eof(float),d_pitchC/sizeof(float));
    cudaMalloc部分改成：
      float* d_A;
cutilSafeCall(cudaMallocPitch((void**)&d_A,&d_pitchA,
wA*sizeof(float),wB));
      float* d_B;
cutilSafeCall(cudaMallocPitch((void**)&d_B,&d_pitchB,
wB*sizeof(float),wA));
      float* d_C;
cutilSafeCall(cudaMallocPitch((void**)&d_C,&d_pitchC,
wB*sizeof(float),wB));
    矩陣內(nèi)存與顯存之間的讀取都需要做相應的修改：    
cutilSafeCall(cudaMemcpy2D(d_A,d_pitchA,A,wA*sizeof(float),
wA*siz
eof(float),wB,cudaMemcpyHostToDevice));    
cutilSafeCall(cudaMemcpy2D(d_B,d_pitchB,B,wB*sizeof(float),
wB*sizeof(float),wA,cudaMemcpyHostToDevice));
cutilSafeCall(cudaMemcpy2D(C,wB*sizeof(float),d_C,d_pitchC,
wB*sizeof(float),wB,cudaMemcpyDeviceToHost));
    在數(shù)值分析，Kahan求和算法(也稱作補償總和)能顯著減少浮點數(shù)運算的誤差，在CUDA矩陣乘法中可以通過使用Kahan求和算法來提高計算精準度[5]。算法如下：
    for (int k = 0; k < BLOCK_DIM; ++k)
    {
            float t;
        comp -= AS[ty][k] * BS[k][tx];
        t = Csub - comp;
        comp = (t - Csub) + comp;
        Csub = t;
    }
3 測試環(huán)境及實驗結(jié)果
    測試的硬件環(huán)境：CPU使用的是AMD Athlon II X2 245處理器,核心數(shù)為2，該處理器主頻為2.9 GHz，峰值運算能力約為17.4 GFLOPS；GPU使用的是NVIDIA GeForce 9800M GTS，有8個SM即有64個SP單元，顯存帶寬為51.2 GB/s，GPU核心頻率為0.625 GHz，單精度浮點計算能力為240 GFLOPS，屬于NVIDIA中端顯卡。測試的軟件環(huán)境：Windows XP系統(tǒng)，CUDA toolkit 3.0,Visual Studio 2008，CUDA計算能力為1.1。
    在程序運行的測試中，對矩陣規(guī)模由256×256～2 048×2 048逐漸增大，實驗數(shù)據(jù)均是三次測試取得的平均值，這樣實驗的結(jié)果更準確。加速比是指程序在CPU上運行的時間與程序在GPU上運行所需的時間之比。峰值比是指運算速度與GPU單精度浮點運算能力之比。最后求得在各種矩陣規(guī)模運行下的加速比及峰值比。實驗結(jié)果如表1所示。

    實驗結(jié)果表明：當矩陣維數(shù)小于320×320時，帶狀劃分加速比小于1，說明CPU運算時間要小于一維帶狀劃分時GPU的運算時間，這說明GPU計算時，從內(nèi)存復制矩陣到顯存和把結(jié)果矩陣從顯存拷貝回內(nèi)存過程中消耗了一些時間[6]。隨著矩陣維數(shù)的增大，CPU的運算時間呈現(xiàn)級數(shù)增長，而GPU運算時間只是小幅度增長。此時GPU強大的浮點運算能力凸顯出來，加速比在矩陣維數(shù)為2 048時最大為1 079.64，CPU上Intel MKL矩陣乘法比文中所用的CPU矩陣乘法快了200多倍，但是依靠GPU流多處理的并行執(zhí)行能力，GPU上的實現(xiàn)方法還是比Intel MKL快了5倍左右。運用CUDA的軟硬件架構(gòu)使得GPU合理組織數(shù)據(jù)，使得內(nèi)存的讀取節(jié)省了大量時間。峰值比也有很大的提高，峰值比說明了算法對GPU強大浮點運算能力的利用，對GPU相應算法的對比具有很高的參考價值。
    通過矩陣乘法在CPU與GPU上不同的性能表現(xiàn)可以發(fā)現(xiàn)，NVIDIA公司推出的CUDA使某些大運算量的計算可以從大型計算機或者超級計算機轉(zhuǎn)移到個人PC，這一新技術不僅使科研縮減了成本，同時也為科學領域進行大規(guī)模運算提供了新方法[7]。對于它的未來值得期待，畢竟CUDA已經(jīng)在影視制作、計算金融、流體力學、醫(yī)學成像、石油天然氣數(shù)據(jù)收集、地質(zhì)勘探及超級計算機的建立等領域取得了成功。
參考文獻
[1] NVIDIA Corporation.NVIDIA CUDA Programming Guide  Version3.0[EB/OL].(2010-02-10)[2011-08-20].http://cuda.csdn.net/.
[2] 張舒，褚艷利，趙開勇，等.GPU高性能并行運算之CUDA[M].北京：中國水利水電出版社，2009.
[3] Ye Zhenyu.GPU assignment 5KK70[DB/OL].(2009-11-05)[2011-09-01].http://wenku.baidu.com/view/9cd2e372027-68e9951e738e5.html.
[4] NVIDIA Corporation.NVIDIA CUDA CUBLAS library PG-00000-002_V3.0[EB/OL].(2010-02-10)[2011-09-10].http://cuda.csdn.net/.
[5] Hotball.深入淺出談CUDA技術[DB/OL].(2008-11-21) [2011-09-15].http://www.pcinlife.com/article/graphics/2008-06-04/1212575164d532_3.html.
[6] 劉進鋒，郭雷.CPU與GPU上幾種矩陣乘法的比較與分析[J].計算機工程與應用，2011,47(19)：9-23.
[7] 肖江，胡柯良，鄧元勇.基于CUDA的矩陣乘法和FFT性能測試[J].計算機工程，2009.35（10）：7-10.