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

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

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

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

    因?yàn)槠灞P劃分中涉及到的是二維數(shù)組，cudaMalloc2D()能確保分配二維數(shù)組并且能分配適當(dāng)?shù)奶畛湟詽M足對(duì)齊要求，還能確保在訪問(wèn)行地址或者二維數(shù)組與其他設(shè)備內(nèi)存之間的數(shù)據(jù)復(fù)制能達(dá)到最佳性能。
    二維棋盤劃分方法僅限于數(shù)組大小必須是BLOCK_SIZE的整數(shù)倍，若矩陣維數(shù)并不是16的整數(shù)倍，則會(huì)造成運(yùn)算效率的下降，此時(shí)可以利用CUDA架構(gòu)特點(diǎn)和CUDA提供的cudaMallocPitch()函數(shù)來(lái)解決此問(wèn)題。cudaMallocPitch()可以自動(dòng)地以最佳倍數(shù)來(lái)分配內(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)存與顯存之間的讀取都需要做相應(yīng)的修改：    
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求和算法(也稱作補(bǔ)償總和)能顯著減少浮點(diǎn)數(shù)運(yùn)算的誤差，在CUDA矩陣乘法中可以通過(guò)使用Kahan求和算法來(lái)提高計(jì)算精準(zhǔn)度[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 測(cè)試環(huán)境及實(shí)驗(yàn)結(jié)果
    測(cè)試的硬件環(huán)境：CPU使用的是AMD Athlon II X2 245處理器,核心數(shù)為2，該處理器主頻為2.9 GHz，峰值運(yùn)算能力約為17.4 GFLOPS；GPU使用的是NVIDIA GeForce 9800M GTS，有8個(gè)SM即有64個(gè)SP單元，顯存帶寬為51.2 GB/s，GPU核心頻率為0.625 GHz，單精度浮點(diǎn)計(jì)算能力為240 GFLOPS，屬于NVIDIA中端顯卡。測(cè)試的軟件環(huán)境：Windows XP系統(tǒng)，CUDA toolkit 3.0,Visual Studio 2008，CUDA計(jì)算能力為1.1。
    在程序運(yùn)行的測(cè)試中，對(duì)矩陣規(guī)模由256×256～2 048×2 048逐漸增大，實(shí)驗(yàn)數(shù)據(jù)均是三次測(cè)試取得的平均值，這樣實(shí)驗(yàn)的結(jié)果更準(zhǔn)確。加速比是指程序在CPU上運(yùn)行的時(shí)間與程序在GPU上運(yùn)行所需的時(shí)間之比。峰值比是指運(yùn)算速度與GPU單精度浮點(diǎn)運(yùn)算能力之比。最后求得在各種矩陣規(guī)模運(yùn)行下的加速比及峰值比。實(shí)驗(yàn)結(jié)果如表1所示。

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