當(dāng)前位置：
首頁(yè) >
深入浅出谈CUDA（二）
發(fā)布時(shí)間：2025/3/17
41
豆豆
生活随笔
收集整理的這篇文章主要介紹了
深入浅出谈CUDA（二）
小編覺(jué)得挺不錯(cuò)的,現(xiàn)在分享給大家,幫大家做個(gè)參考.
前面介紹的計(jì)算平方和的程序，似乎沒(méi)有什么實(shí)用價(jià)值。所以我們的第二個(gè)?CUDA?程序，要做一個(gè)確實(shí)有（某些）實(shí)用價(jià)值的程序，也就是進(jìn)行矩陣乘法。而且，這次我們會(huì)使用浮點(diǎn)數(shù)。
雖然矩陣乘法有點(diǎn)老套，不過(guò)因?yàn)樗喈?dāng)簡(jiǎn)單，而且也可以用來(lái)介紹一些有關(guān)?CUDA?的有趣性質(zhì)。
為了單純起見(jiàn)，我們這里以方形的矩陣為例子。基本上，假設(shè)有兩個(gè)矩陣?A?和?B，則計(jì)算?AB?=?C?的方法如下：????for(i?=?0;?i?<?n;?i++)?{
????????for(j?=?0;?j?<?n;?j++)?{
????????????C[i][j]?=?0;
????????????for(k?=?0;?k?<?n;?k++)?{
????????????????C[i][j]?+=?A[i][k]?*?B[k][j];
????????????}
????????}
????}
一開(kāi)始，我們先準(zhǔn)備好產(chǎn)生數(shù)據(jù)、設(shè)定?CUDA?等等的工作：????int?main()
????{
????????float?*a,?*b,?*c,?*d;
????????int?n?=?1000;
????????if(!InitCUDA())?return?0;
????????a?=?(float*)?malloc(sizeof(float)?*?n?*?n);
????????b?=?(float*)?malloc(sizeof(float)?*?n?*?n);
????????c?=?(float*)?malloc(sizeof(float)?*?n?*?n);
????????d?=?(float*)?malloc(sizeof(float)?*?n?*?n);
????????srand(0);
????????matgen(a,?n,?n);
????????matgen(b,?n,?n);
????????clock_t?time?=?matmultCUDA(a,?n,?b,?n,?c,?n,?n);
????????matmult(a,?n,?b,?n,?d,?n,?n);
????????compare_mat(c,?n,?d,?n,?n);
????????double?sec?=?(double)?time?/?CLOCKS_PER_SEC;
????????printf("Time?used:?%.2f?(%.2lf?GFLOPS)\n",?sec,
????????????2.0?*?n?*?n?*?n?/?(sec?*?1E9));
????????return?0;
????}InitCUDA?函式和第一個(gè)?CUDA?程序一樣，可以直接參考前面的文章。以下是上面用到的一些其它的函式：產(chǎn)生矩陣：????void?matgen(float*?a,?int?lda,?int?n)
????{
????????int?i,?j;
????????for(i?=?0;?i?<?n;?i++)?{
????????????for(j?=?0;?j?<?n;?j++)?{
????????????????a[i?*?lda?+?j]?=?(float)?rand()?/?RAND_MAX?+
?????????????????????(float)?rand()?/?(RAND_MAX?*?RAND_MAX);
????????????}
????????}
????}這個(gè)函式只是利用隨機(jī)數(shù)生成器把矩陣填滿?0?~?1?之間的數(shù)字。特別注意到因?yàn)?C?語(yǔ)言中無(wú)法聲明變動(dòng)大小的二維矩陣，所以我們使用?i?*?lda?+?j?的方式。進(jìn)行矩陣乘法：????void?matmult(const?float*?a,?int?lda,?const?float*?b,?int?ldb,
????????float*?c,?int?ldc,?int?n)
????{
????????int?i,?j,?k;
????????for(i?=?0;?i?<?n;?i++)?{
????????????for(j?=?0;?j?<?n;?j++)?{
????????????????double?t?=?0;
????????????????for(k?=?0;?k?<?n;?k++)?{
????????????????????t?+=?a[i?*?lda?+?k]?*?b[k?*?ldb?+?j];
????????????????}
????????????????c[i?*?ldc?+?j]?=?t;
????????????}
????????}
????}這是以?CPU?進(jìn)行矩陣乘法、用來(lái)進(jìn)行驗(yàn)證答案正確與否的程序。特別注意到它用?double?來(lái)儲(chǔ)存暫時(shí)的計(jì)算結(jié)果，以提高精確度。
驗(yàn)證結(jié)果：????void?compare_mat(const?float*?a,?int?lda,
????????const?float*?b,?int?ldb,?int?n)
????{
????????float?max_err?=?0;
????????float?average_err?=?0;
????????int?i,?j;
????????for(i?=?0;?i?<?n;?i++)?{
????????????for(j?=?0;?j?<?n;?j++)?{
????????????????if(b[i?*?ldb?+?j]?!=?0)?{
????????????????????float?err?=?fabs((a[i?*?lda?+?j]?-
?????????????????????????b[i?*?ldb?+?j])?/?b[i?*?ldb?+?j]);
????????????????????if(max_err?<?err)?max_err?=?err;
????????????????????average_err?+=?err;
????????????????}
????????????}
????????}
????????printf("Max?error:?%g?Average?error:?%g\n",
????????????max_err,?average_err?/?(n?*?n));
????}這個(gè)函式計(jì)算兩個(gè)矩陣的最大相對(duì)誤差和平均相對(duì)誤差，并把結(jié)果印出來(lái)。最后是?CUDA?的矩陣乘法的部份：????#define?NUM_THREADS?256
????clock_t?matmultCUDA(const?float*?a,?int?lda,
????????const?float*?b,?int?ldb,?float*?c,?int?ldc,?int?n)
????{
????????float?*ac,?*bc,?*cc;
????????clock_t?start,?end;
????????start?=?clock();
????????cudaMalloc((void**)?&ac,?sizeof(float)?*?n?*?n);
?????????cudaMalloc((void**)?&bc,?sizeof(float)?*?n?*?n);
????????cudaMalloc((void**)?&cc,?sizeof(float)?*?n?*?n);
????????cudaMemcpy2D(ac,?sizeof(float)?*?n,?a,?sizeof(float)?*?lda,
????????????sizeof(float)?*?n,?n,?cudaMemcpyHostToDevice);
????????cudaMemcpy2D(bc,?sizeof(float)?*?n,?b,?sizeof(float)?*?ldb,
????????????sizeof(float)?*?n,?n,?cudaMemcpyHostToDevice);
????????int?blocks?=?(n?+?NUM_THREADS?-?1)?/?NUM_THREADS;
????????matMultCUDA<<<blocks?*?n,?NUM_THREADS>>>
????????????(ac,?n,?bc,?n,?cc,?n,?n);
????????cudaMemcpy2D(c,?sizeof(float)?*?ldc,?cc,?sizeof(float)?*?n,
????????sizeof(float)?*?n,?n,?cudaMemcpyDeviceToHost);
????????cudaFree(ac);
????????cudaFree(bc);
????????cudaFree(cc);
????????end?=?clock();
????????return?end?-?start;
????}這個(gè)函式相當(dāng)單純，就是在顯卡內(nèi)存中配置存放矩陣的內(nèi)存，然后把主內(nèi)存中的矩陣數(shù)據(jù)復(fù)制到顯卡內(nèi)存上。不過(guò)，因?yàn)槲覀兊木仃嚦朔ê娇梢灾付?pitch（即?lda、ldb、和?ldc），所以如果用一般的?cudaMemcpy?函式來(lái)復(fù)制內(nèi)存的話，會(huì)需要每個(gè)?row?都分開(kāi)復(fù)制，那會(huì)需要呼叫很多次?cudaMemcpy?函式，會(huì)使效率變得很差。因此，在這里我們用了一個(gè)新的?cudaMemcpy2D?函式，它是用來(lái)復(fù)制二維數(shù)組，可以指定數(shù)組的?pitch。這樣就可以透過(guò)一次函數(shù)調(diào)用就可以了。進(jìn)行計(jì)算的?kernel?如下：????__global__?static?void?matMultCUDA(const?float*?a,?size_t?lda,
????????const?float*?b,?size_t?ldb,?float*?c,?size_t?ldc,?int?n)
????{
????????const?int?tid?=?threadIdx.x;
????????const?int?bid?=?blockIdx.x;
????????const?int?idx?=?bid?*?blockDim.x?+?tid;
????????const?int?row?=?idx?/?n;
????????const?int?column?=?idx?%?n;
????????int?i;
????????if(row?<?n?&&?column?<?n)?{
????????????float?t?=?0;
????????????for(i?=?0;?i?<?n;?i++)?{
????????????????t?+=?a[row?*?lda?+?i]?*?b[i?*?ldb?+?column];
????????????}
????????????c[row?*?ldc?+?column]?=?t;
????????}
????}這個(gè)函式一開(kāi)始先從?bid?和?tid?計(jì)算出這個(gè)?thread?應(yīng)該計(jì)算的?row?和?column，在判斷?row?和?column?在范圍內(nèi)之后，就直接進(jìn)行計(jì)算，并把結(jié)果寫到?c?矩陣中，是非常單純的函式。在?GeForce?8800GT?上實(shí)際執(zhí)行的結(jié)果如下：????Max?error:?2.01484e-006?Average?error:?3.36637e-007
????Time?used:?1.1560?(1.73?GFLOPS)可以看到兩個(gè)問(wèn)題：1?很明顯的，執(zhí)行效率相當(dāng)?shù)吐洹?/span>2?最大相對(duì)誤差偏高。理想上應(yīng)該要低于?1e-6。計(jì)算結(jié)果的誤差偏高的原因是，在?CPU?上進(jìn)行計(jì)算時(shí)，我們使用?double（即?64?bits?浮點(diǎn)數(shù)）來(lái)累進(jìn)計(jì)算過(guò)程，而在?GPU?上則只能用?float（32?bits?浮點(diǎn)數(shù)）。在累加大量數(shù)字的時(shí)候，由于累加結(jié)果很快會(huì)變大，因此后面的數(shù)字很容易被舍去過(guò)多的位數(shù)。由于?CUDA?的浮點(diǎn)數(shù)運(yùn)算，在進(jìn)行加、減、乘法時(shí)是符合?IEEE?754?規(guī)定的精確度的，因此，我們可以利用?Kahan's?Summation?Formula?來(lái)提高精確度。把程序改成：????if(row?<?n?&&?column?<?n)?{
????????float?t?=?0;
????????float?y?=?0;
????????for(i?=?0;?i?<?n;?i++)?{
????????????float?r;
????????????y?-=?a[row?*?lda?+?i]?*?b[i?*?ldb?+?column];
????????????r?=?t?-?y;
????????????y?=?(r?-?t)?+?y;
????????????t?=?r;
????????}
????}修改后的程序的執(zhí)行結(jié)果是：????Max?error:?1.19209e-007?Average?error:?4.22751e-008
????Time?used:?1.1560?(1.73?GFLOPS)可以看到相對(duì)誤差有很大的改善，效率則沒(méi)什么變化。由于?Kahan's?Summation?Formula?需要的運(yùn)算量提高，但是效率卻沒(méi)有什么改變，可以看出這個(gè)?kernel?主要的瓶頸應(yīng)該是在內(nèi)存的存取動(dòng)作上。這是因?yàn)橛写罅康膬?nèi)存讀取是重復(fù)的。例如，矩陣?a?的一個(gè)?row?在每次進(jìn)行計(jì)算時(shí)都被重復(fù)讀入，但這是相當(dāng)浪費(fèi)的。這樣的計(jì)算方式，總共需要讀取?2*n3?次內(nèi)存。如果讓一個(gè)?row?只需要讀入一次的話，就可以減到為?n3+n2?次。
?和我們的第一個(gè)?CUDA?程序一樣，我們可以利用?shared?memory?來(lái)儲(chǔ)存每個(gè)?row?的數(shù)據(jù)。不過(guò)，因?yàn)橹挥型粋€(gè)?block?的?threads?可以共享?shared?memory，因此現(xiàn)在一個(gè)?row?只能由同一個(gè)?block?的?threads?來(lái)進(jìn)行計(jì)算。另外我們也需要能存放一整個(gè)?row?的?shared?memory。因此，把先把呼叫?kernel?的部份改成：????????matMultCUDA<<<n,?NUM_THREADS,?sizeof(float)?*?n>>>
????????????(ac,?n,?bc,?n,?cc,?n,?n);kernel?的部份則改成：????__global__?static?void?matMultCUDA(const?float*?a,?size_t?lda,
????????const?float*?b,?size_t?ldb,?float*?c,?size_t?ldc,?int?n)
????{
????????extern?__shared__?float?data[];
????????const?int?tid?=?threadIdx.x;
????????const?int?row?=?blockIdx.x;
????????int?i,?j;
????????for(i?=?tid;?i?<?n;?i?+=?blockDim.x)?{
????????????data[i]?=?a[row?*?lda?+?i];
????????}
????????__syncthreads();
????????for(j?=?tid;?j?<?n;?j?+=?blockDim.x)?{
????????????float?t?=?0;
????????????float?y?=?0;
????????????for(i?=?0;?i?<?n;?i++)?{
????????????????float?r;
????????????????y?-=?data[i]?*?b[i?*?ldb?+?j];
????????????????r?=?t?-?y;
????????????????y?=?(r?-?t)?+?y;
????????????????t?=?r;
????????????}
????????????c[row?*?ldc?+?j]?=?t;
????????}
????}
第一個(gè)部份先把整個(gè)?row?讀到?shared?memory?中，而第二個(gè)部份則進(jìn)行計(jì)算，并沒(méi)有太大的變化。主要的差別是現(xiàn)在一個(gè)?row?只由一個(gè)?block?進(jìn)行計(jì)算。在?GeForce?8800GT?上，執(zhí)行的結(jié)果是：????Max?error:?1.19209e-007??Average?error:?4.22751e-008
????Time?used:?0.4220???(4.74?GFLOPS)很明顯的，計(jì)算的結(jié)果并沒(méi)有改變，不過(guò)速度則提高了超過(guò)一倍。雖然如此，但是這樣的效率仍不盡理想，因?yàn)槔碚撋?GeForce?8800GT?有超過(guò)?300GFLOPS?的運(yùn)算性能。即使是把?Kahan's?Summation?Formula?所需要的額外運(yùn)算考慮進(jìn)去，這樣的效率仍然連理論最大值的十分之一都不到。會(huì)有這樣的結(jié)果，原因其實(shí)還是同樣的：對(duì)內(nèi)存的存取次數(shù)太多了。雖然現(xiàn)在?A?矩陣的?row?的數(shù)據(jù)已經(jīng)不再需要重復(fù)讀取，但是?B?矩陣的?column?的數(shù)據(jù)仍然一直被重復(fù)讀取。另一個(gè)問(wèn)題比較不是那么明顯：對(duì)?B?矩陣的讀取，雖然看起來(lái)不連續(xù)，但實(shí)際上它是連續(xù)的。這是因?yàn)椴煌?thread?會(huì)讀取不同的?column，因此同時(shí)間每個(gè)?thread?讀取的各個(gè)?column?加起來(lái)，就是一個(gè)連續(xù)的內(nèi)存區(qū)塊。那么，為什么效率還是不佳呢？這是因?yàn)?#xff0c;GPU?上的內(nèi)存控制器，從某個(gè)固定的倍數(shù)地址開(kāi)始讀取，才會(huì)有最高的效率（例如?16?bytes?的倍數(shù)）。由于矩陣大小并不是?16?的倍數(shù)（這里使用的是?1000x1000?的矩陣），所以造成效率不佳的情形。要解決這個(gè)問(wèn)題，我們可以在?cudaMalloc的時(shí)候稍微修改一下，讓寬度變成?適當(dāng)?shù)谋稊?shù)就可以了。但是，適當(dāng)?shù)谋稊?shù)是多少呢？幸運(yùn)的是，我們并不需要知道這些細(xì)節(jié)。CUDA?提供了一個(gè)?cudaMallocPitch的函式，可以自動(dòng)以最佳的倍數(shù)來(lái)配置內(nèi)存。因此，我們可以把?cudaMalloc?的部份改成：????size_t?pitch_a,?pitch_b,?pitch_c;
????cudaMallocPitch((void**)?&ac,?&pitch_a,?sizeof(float)?*?n,?n);
????cudaMallocPitch((void**)?&bc,?&pitch_b,?sizeof(float)?*?n,?n);
????cudaMallocPitch((void**)?&cc,?&pitch_c,?sizeof(float)?*?n,?n);cudaMallocPitch?函式會(huì)以適當(dāng)?shù)谋稊?shù)配置內(nèi)存，并把配置的寬度傳回。因此，在把矩陣復(fù)制到顯卡內(nèi)存上時(shí)，要使用它傳回的寬度：????cudaMemcpy2D(ac,?pitch_a,?a,?sizeof(float)?*?lda,
????????sizeof(float)?*?n,?n,?cudaMemcpyHostToDevice);
????cudaMemcpy2D(bc,?pitch_b,?b,?sizeof(float)?*?ldb,
????????sizeof(float)?*?n,?n,?cudaMemcpyHostToDevice);呼叫?kernel?的部份也需要修改：????matMultCUDA<<<n,?NUM_THREADS,?sizeof(float)?*?n>>>
????????(ac,?pitch_a?/?sizeof(float),?bc,?pitch_b?/?sizeof(float),
????????cc,?pitch_c?/?sizeof(float),?n);同樣的，把計(jì)算結(jié)果復(fù)制回到主內(nèi)存時(shí)，也要使用傳回的寬度值：????cudaMemcpy2D(c,?sizeof(float)?*?ldc,?cc,?pitch_c,
????????sizeof(float)?*?n,?n,?cudaMemcpyDeviceToHost);這樣就完成了。Kernel?部份則不需要修改。這樣的修改有多大的效果呢？在?GeForce?8800GT?上執(zhí)行，結(jié)果如下：????Max?error:?1.19209e-007??Average?error:?4.22751e-008
????Time?used:?0.1250???(16.00?GFLOPS)可以看到，執(zhí)行速度又再大幅提高了三倍多！而這只是把內(nèi)存的配置方式稍微修改一下而已。雖然執(zhí)行速度提高了很多，但是，和前面提到的理論值相比，其實(shí)還是有相當(dāng)?shù)牟罹?。這是因?yàn)?#xff0c;前面也提到過(guò)，這樣的做法需要?n3+n2次的內(nèi)存讀取，和?n2次?的內(nèi)存寫入動(dòng)作。由于?n?=?1000，每個(gè)數(shù)字的大小是?32?bits，所以總共的內(nèi)存存取數(shù)據(jù)量約為?4GB。除以實(shí)際執(zhí)行的時(shí)間?0.125?秒，得到的帶寬數(shù)值是約?32GB/s，這已經(jīng)接近?GeForce?8800GT?顯卡內(nèi)存的帶寬了。由于我們計(jì)算時(shí)間的時(shí)候，把配置內(nèi)存、以及數(shù)據(jù)的復(fù)制動(dòng)作也計(jì)算進(jìn)去，因此實(shí)際上花費(fèi)在?kernel?的時(shí)間是更短的（約?0.09?秒）。因此，可以很明顯的看出，這個(gè)程序的效率，是受限于內(nèi)存帶寬的。
上一節(jié)的結(jié)論顯示出，矩陣乘法的程序，效率是受限于內(nèi)存帶寬的。那有沒(méi)有辦法降低內(nèi)存的存取次數(shù)呢？答案當(dāng)然是有的，不然就不會(huì)有這一節(jié)了?:)要進(jìn)一步降低內(nèi)存帶寬的使用，可以注意到，在上一節(jié)的方法中，雖然?A?矩陣的存取次數(shù)被減至最低，但是?B?矩陣的存取次數(shù)并沒(méi)有減少。這是因?yàn)槲覀冎粚?A?矩陣的?row?加載到?shared?memory?中，但是?B?矩陣的?column?也是有被重復(fù)使用的。理想上應(yīng)該也可以避免重復(fù)加載才對(duì)。不過(guò)，由于?B?矩陣的?column?使用的時(shí)機(jī)，和?A?矩陣的?row?是不同的，所以并不能直接這樣做。解決方法是?"blocking"。也就是把整個(gè)矩陣乘法的動(dòng)作，切割成很多小矩陣的乘法。例如，要計(jì)算?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)?+?...這樣一來(lái)，我們就可以把兩個(gè)小矩陣加載到?shared?memory，則小矩陣本身的乘法就不需要再存取任何外部的內(nèi)存了！這樣一來(lái)，假設(shè)小矩陣的大小是?k，則實(shí)際上需要的內(nèi)存存取次數(shù)就會(huì)變成約?2k2(n/k)3?=?2n3/k。由于目前?CUDA?每個(gè)?block?的?thread?數(shù)目最多是?512，因此?k?=?16?似乎是一個(gè)相當(dāng)理想的數(shù)字（共?256?個(gè)?threads）。因此，對(duì)于一個(gè)?n?=?1000?的矩陣來(lái)說(shuō)，我們可以把內(nèi)存存取的量減少到約?500MB，也就是上一節(jié)的存取量的?1/8。理論上，這樣應(yīng)該可以讓效率提高八倍（假設(shè)沒(méi)有遇到別的瓶頸）。為了方便進(jìn)行區(qū)塊的計(jì)算，我們讓每個(gè)?block?有?16x16?個(gè)?threads，再建立?(n/16)x(n/16)?個(gè)?blocks。把呼叫?kernel?的地方改成：????int?bx?=?(n?+?BLOCK_SIZE?-?1)?/?BLOCK_SIZE;
????dim3?blocks(bx,?bx);
????dim3?threads(BLOCK_SIZE,?BLOCK_SIZE);
????matMultCUDA<<<blocks,?threads>>>(ac,?pitch_a?/?sizeof(float),
????????bc,?pitch_b?/?sizeof(float),?cc,?pitch_c?/?sizeof(float),?n);BLOCK_SIZE?則是定義成?16。dim3?是?CUDA?的一種數(shù)據(jù)型態(tài)，表示一個(gè)?3D?的向量。在這里，我們透過(guò)?dim3?來(lái)建立?16x16?個(gè)?threads?的?block，和?(n/16)x(n/16)?個(gè)?blocks。Kernel?程序的部份，則改成：????__global__?static?void?matMultCUDA(const?float*?a,?size_t?lda,
????????const?float*?b,?size_t?ldb,?float*?c,?size_t?ldc,?int?n)
????{
????????__shared__?float?matA[BLOCK_SIZE][BLOCK_SIZE];
????????__shared__?float?matB[BLOCK_SIZE][BLOCK_SIZE];
????????const?int?tidc?=?threadIdx.x;
????????const?int?tidr?=?threadIdx.y;
????????const?int?bidc?=?blockIdx.x?*?BLOCK_SIZE;
????????const?int?bidr?=?blockIdx.y?*?BLOCK_SIZE;
????????int?i,?j;
????????float?results?=?0;
????????float?comp?=?0;
????????for(j?=?0;?j?<?n;?j?+=?BLOCK_SIZE)?{
????????????if(tidr?+?bidr?<?n?&&?tidc?+?j?<?n)?{
????????????????matA[tidr][tidc]?=?a[(tidr?+?bidr)?*?lda?+?tidc?+?j];
????????????}
????????????else?{
????????????????matA[tidr][tidc]?=?0;
????????????}
????????????if(tidr?+?j?<?n?&&?tidc?+?bidc?<?n)?{
????????????????matB[tidr][tidc]?=?b[(tidr?+?j)?*?ldb?+?tidc?+?bidc];
????????????}
????????????else?{
????????????????matB[tidr][tidc]?=?0;
????????????}
????????????__syncthreads();
????????????for(i?=?0;?i?<?BLOCK_SIZE;?i++)?{
????????????????float?t;
????????????????comp?-=?matA[tidr][i]?*?matB[i][tidc];
????????????????t?=?results?-?comp;
????????????????comp?=?(t?-?results)?+?comp;
????????????????results?=?t;
????????????}
????????????__syncthreads();
????????}
????????if(tidr?+?bidr?<?n?&&?tidc?+?bidc?<?n)?{
????????????c[(tidr?+?bidr)?*?ldc?+?tidc?+?bidc]?=?results;
????????}
????}注意到因?yàn)槲覀儸F(xiàn)在使用?16x16?的?threads，因此?threadIdx?變量可以取得?threadIdx.x?和?threadIdx.y，范圍分別是?0?~?15。blockIdx.x?和?blockIdx.y?變量也是同樣的情形，范圍分別是?0?~?n/16。在程序中，因?yàn)榫仃嚨拇笮〔灰欢〞?huì)是?16?的倍數(shù)，因此需要使用?if?判斷式檢查是否超出矩陣范圍。這個(gè)版本在?GeForce?8800GT?上的執(zhí)行結(jié)果如下：????Max?error:?1.19209e-007??Average?error:?4.22751e-008
????Time?used:?0.0780???(25.64?GFLOPS)速度雖然提高了，但是似乎并沒(méi)有達(dá)到預(yù)期中的八倍。當(dāng)然，前面提到過(guò)，我們?cè)谟?jì)算時(shí)間時(shí)，把一些復(fù)制內(nèi)存、配置內(nèi)存的動(dòng)作也計(jì)算在內(nèi)，這些動(dòng)作的時(shí)?間并不會(huì)縮短。實(shí)際上?kernel?的運(yùn)行時(shí)間，大約是?0.053?秒左右（約略相當(dāng)于?38GFLOPS），比上一節(jié)的版本快了將近一倍。如果這一版程序已經(jīng)不再限于內(nèi)存帶寬，那為什么沒(méi)有達(dá)到預(yù)期的效率呢？這是因?yàn)檫@一版程序已經(jīng)是限于指令周期了。除了使用?Kahan's?Summation?Formula?會(huì)需要更多的運(yùn)算之外，程序中也有大量計(jì)算矩陣地址的乘法等等，這都會(huì)需要花費(fèi)運(yùn)算資源。另外，那些用來(lái)判斷超出矩陣范圍的?if?判斷式，也會(huì)有一定的影響。要把那些?if?判斷式去掉，有一個(gè)方法是，在配置內(nèi)存時(shí)，就配置成?16?的倍數(shù)，并在復(fù)制矩陣到顯卡內(nèi)存之前，先將它清為?0。如下所示：????int?newn?=?((n?+?BLOCK_SIZE?-?1)?/?BLOCK_SIZE)?*?BLOCK_SIZE;
????cudaMallocPitch((void**)?&ac,?&pitch_a,
?????????sizeof(float)?*?newn,?newn);
????cudaMallocPitch((void**)?&bc,?&pitch_b,
????????sizeof(float)?*?newn,?newn);
????cudaMallocPitch((void**)?&cc,?&pitch_c,
????????sizeof(float)?*?newn,?newn);
????cudaMemset(ac,?0,?pitch_a?*?newn);
????cudaMemset(bc,?0,?pitch_b?*?newn);?這樣一來(lái)，我們就可以把?kernel?中的?if?判斷式都移除了：????__global__?static?void?matMultCUDA(const?float*?a,?size_t?lda,
????????const?float*?b,?size_t?ldb,?float*?c,?size_t?ldc,?int?n)
????{
????????__shared__?float?matA[BLOCK_SIZE][BLOCK_SIZE];
????????__shared__?float?matB[BLOCK_SIZE][BLOCK_SIZE];
????????const?int?tidc?=?threadIdx.x;
????????const?int?tidr?=?threadIdx.y;
????????const?int?bidc?=?blockIdx.x?*?BLOCK_SIZE;
????????const?int?bidr?=?blockIdx.y?*?BLOCK_SIZE;
????????int?i,?j;
????????float?results?=?0;
????????float?comp?=?0;
????????for(j?=?0;?j?<?n;?j?+=?BLOCK_SIZE)?{
????????????matA[tidr][tidc]?=?a[(tidr?+?bidr)?*?lda?+?tidc?+?j];
????????????matB[tidr][tidc]?=?b[(tidr?+?j)?*?ldb?+?tidc?+?bidc];
????????????__syncthreads();
????????????for(i?=?0;?i?<?BLOCK_SIZE;?i++)?{
????????????????float?t;
????????????????comp?-=?matA[tidr][i]?*?matB[i][tidc];
????????????????t?=?results?-?comp;
????????????????comp?=?(t?-?results)?+?comp;
????????????????results?=?t;
????????????}
????????????__syncthreads();
????????}
????????c[(tidr?+?bidr)?*?ldc?+?tidc?+?bidc]?=?results;
????}這個(gè)版本的執(zhí)行結(jié)果是：????Max?error:?1.19209e-007??Average?error:?4.22751e-008
????Time?used:?0.0780???(25.64?GFLOPS)似乎沒(méi)有改善。不過(guò)，實(shí)際上?kernel?的運(yùn)行時(shí)間已經(jīng)減少到?0.042?秒（約略相當(dāng)于?48GFLOPS）。
有些讀者可能會(huì)想，如果把?block?再變得更大（例如?32x32）是否會(huì)有幫助呢？當(dāng)然，由于最后的程序已經(jīng)不再是受限于內(nèi)存帶寬（在?0.042?秒內(nèi)存取?500MB?的數(shù)據(jù)約相當(dāng)于?12GB/s?的帶寬），所以把?block?再加大并不會(huì)有幫助了。而且，由于一個(gè)?block?內(nèi)的?thread?數(shù)目最多只能到?512?個(gè)，將?block?變大也會(huì)造成很多額外負(fù)擔(dān)。而且?shared?memory?的大小也有限制（GeForce?8800GT?的?shared?memory?大小限制是?16384?bytes），所以也不能任意增加?block?的大小。最后一版程序的完整檔案可以從這里下載。GPU?的硬件架構(gòu)?這里我們會(huì)簡(jiǎn)單介紹，NVIDIA?目前支持?CUDA?的?GPU，其在執(zhí)行?CUDA?程序的部份（基本上就是其?shader?單元）的架構(gòu)。這里的數(shù)據(jù)是綜合?NVIDIA?所公布的信息，以及?NVIDIA?在各個(gè)研討會(huì)、學(xué)校課程等所提供的數(shù)據(jù)，因此有可能會(huì)有不正確的地方。主要的數(shù)據(jù)源包括?NVIDIA?的?CUDA?Programming?Guide?1.1、NVIDIA?在?Supercomputing?'07?介紹?CUDA?的?session，以及?UIUC?的?CUDA?課程。
目前?NVIDIA?推出的顯示芯片，支持?CUDA?的是?G80?系列的顯示芯片。其中?G80?顯示芯片支持?CUDA?1.0?版，而?G84、G86、G92、G94、G96?則支援?CUDA?1.1?版?；旧?#xff0c;除了最早的?GeForce?8800?Ultra/GTX?及?320MB/640MB?版本的?GeForce?8800GTS、Tesla?等顯卡是?CUDA?1.0?版之外，其它?GeForce?8?系列及?9?系列顯卡都支持?CUDA?1.1。詳細(xì)情形可以參考?CUDA?Programming?Guide?1.1?的?Appendix?A。所有目前支持?CUDA?的?NVIDIA?顯示芯片，其?shader?部份都是由多個(gè)?multiprocessors?組成。每個(gè)?multiprocessor?里包含了八個(gè)?stream?processors，?其組成是四個(gè)四個(gè)一組，也就是說(shuō)實(shí)際上可以看成是有兩組?4D?的?SIMD?處理器。此外，每個(gè)?multiprocessor?還具有?8192?個(gè)寄存器，16KB?的?share?memory，以及?texture?cache?和?constant?cache。大致上如下圖所示：
詳細(xì)的?multiprocessor?信息，都可以透過(guò)?CUDA?的?cudaGetDeviceProperties()?函式或?cuDeviceGetProperties()?函式取得。不過(guò)，目前還沒(méi)有辦法直接取得一個(gè)顯示芯片中有多少?multiprocessor?的信息。
在?CUDA?中，大部份基本的運(yùn)算動(dòng)作，都可以由?stream?processor?進(jìn)行。每個(gè)?stream?processor?都包含一個(gè)?FMA（fused-multiply-add）單元，可以進(jìn)行一個(gè)乘法和一個(gè)加法。比較復(fù)雜的運(yùn)算則會(huì)需要比較長(zhǎng)的時(shí)間。
在執(zhí)行?CUDA?程序的時(shí)候，每個(gè)?stream?processor?就是對(duì)應(yīng)一個(gè)?thread。每個(gè)?multiprocessor?則對(duì)應(yīng)一個(gè)?block。從之前的文章中，可以注意到一個(gè)?block?經(jīng)常有很多個(gè)?thread（例如?256?個(gè)），遠(yuǎn)超過(guò)一個(gè)?multiprocessor?所有的?stream?processor?數(shù)目。這又是怎么回事呢？
實(shí)際上，雖然一個(gè)?multiprocessor?只有八個(gè)?stream?processor，但是由于?stream?processor?進(jìn)行各種運(yùn)算都有?latency，更不用提內(nèi)存存取的?latency，因此?CUDA?在執(zhí)行程序的時(shí)候，是以?warp?為單位。目前的?CUDA?裝置，一個(gè)?warp?里面有?32?個(gè)?threads，分成兩組?16?threads?的?half-warp。由于?stream?processor?的運(yùn)算至少有?4?cycles?的?latency，因此對(duì)一個(gè)?4D?的?stream?processors?來(lái)說(shuō)，一次至少執(zhí)行?16?個(gè)?threads（即?half-warp）才能有效隱藏各種運(yùn)算的?latency。
由于?multiprocessor?中并沒(méi)有太多別的內(nèi)存，因此每個(gè)?thread?的狀態(tài)都是直接保存在?multiprocessor?的寄存器中。所以，如果一個(gè)?multiprocessor?同時(shí)有愈多的?thread?要執(zhí)行，就會(huì)需要愈多的寄存器空間。例如，假設(shè)一個(gè)?block?里面有?256?個(gè)?threads，每個(gè)?thread?用到?20?個(gè)寄存器，那么總共就需要?256x20?=?5,120?個(gè)寄存器才能保存每個(gè)?thread?的狀態(tài)。目前?CUDA?裝置中每個(gè)?multiprocessor?有?8,192?個(gè)寄存器，因此，如果每個(gè)?thread?使用到?16?個(gè)寄存器，那就表示一個(gè)?multiprocessor?同時(shí)最多只能維持?512?個(gè)?thread?的執(zhí)行。如果同時(shí)進(jìn)行的?thread?數(shù)目超過(guò)這個(gè)數(shù)字，那么就會(huì)需要把一部份的數(shù)據(jù)儲(chǔ)存在顯卡內(nèi)存中，就會(huì)降低執(zhí)行的效率了。編者注：在NVIDIA?GT200中的Register?File大小增加了一倍，在FP32下可用的register?file為16K，FP64下是8K。
目前?CUDA?裝置中，每個(gè)?multiprocessor?有?16KB?的?shared?memory。Shared?memory?分成?16?個(gè)?bank。如果同時(shí)每個(gè)?thread?是存取不同的?bank，就不會(huì)產(chǎn)生任何問(wèn)題，存取?shared?memory?的速度和存取寄存器相同。不過(guò)，如果同時(shí)有兩個(gè)（或更多個(gè)）?threads?存取同一個(gè)?bank?的數(shù)據(jù)，就會(huì)發(fā)生?bank?conflict，這些?threads?就必須照順序去存取，而無(wú)法同時(shí)存取?shared?memory?了。Shared?memory?是以?4?bytes?為單位分成?banks。因此，假設(shè)以下的數(shù)據(jù)：????__shared__?int?data[128];那么，data[0]?是?bank?0、data[1]?是?bank?1、data[2]?是?bank?2、…、data[15]?是?bank?15，而?data[16]?又回到?bank?0。由于?warp?在執(zhí)行時(shí)是以?half-warp?的方式執(zhí)行，因此分屬于不同的?half?warp?的?threads，不會(huì)造成?bank?conflict。因此，如果程序在存取?shared?memory?的時(shí)候，使用以下的方式：????int?number?=?data[base?+?tid];那就不會(huì)有任何?bank?conflict，可以達(dá)到最高的效率。但是，如果是以下的方式：????int?number?=?data[base?+?4?*?tid];那么，thread?0?和?thread?4?就會(huì)存取到同一個(gè)?bank，thread?1?和?thread?5?也是同樣，這樣就會(huì)造成?bank?conflict。在這個(gè)例子中，一個(gè)?half?warp?的?16?個(gè)?threads?會(huì)有四個(gè)?threads?存取同一個(gè)?bank，因此存取?share?memory?的速度會(huì)變成原來(lái)的?1/4。一個(gè)重要的例外是，當(dāng)多個(gè)?thread?存取到同一個(gè)?shared?memory?的地址時(shí)，shared?memory?可以將這個(gè)地址的?32?bits?數(shù)據(jù)「廣播」到所有讀取的?threads，因此不會(huì)造成?bank?conflict。例如：????int?number?=?data[3];這樣不會(huì)造成?bank?conflict，因?yàn)樗械?thread?都讀取同一個(gè)地址的數(shù)據(jù)。很多時(shí)候?shared?memory?的?bank?conflict?可以透過(guò)修改數(shù)據(jù)存放的方式來(lái)解決。例如，以下的程序：????data[tid]?=?global_data[tid];
????...
????int?number?=?data[16?*?tid];
會(huì)造成嚴(yán)重的?bank?conflict，為了避免這個(gè)問(wèn)題，可以把數(shù)據(jù)的排列方式稍加修改，把存取方式改成：????int?row?=?tid?/?16;
????int?column?=?tid?%?16;
????data[row?*?17?+?column]?=?global_data[tid];
????...
????int?number?=?data[17?*?tid];這樣就不會(huì)造成?bank?conflict?了。編者注：share?memory在NVIDIA的文檔中其實(shí)還有不同的叫法，例如PDC（Parallel?Data?Cache）、PBSM（per-block?share?memory）。
由于?multiprocessor?并沒(méi)有對(duì)?global?memory?做?cache（如果每個(gè)?multiprocessor?都有自己的?global?memory?cache，將會(huì)需要?cache?coherence?protocol，會(huì)大幅增加?cache?的復(fù)雜度），所以?global?memory?存取的?latency?非常的長(zhǎng)。除此之外，前面的文章中也提到過(guò)?global?memory?的存取，要盡可能的連續(xù)。這是因?yàn)?DRAM?存取的特性所造成的結(jié)果。更精確的說(shuō)，global?memory?的存取，需要是?"coalesced"。所謂的?coalesced，是表示除了連續(xù)之外，而且它開(kāi)始的地址，必須是每個(gè)?thread?所存取的大小的?16?倍。例如，如果每個(gè)?thread?都讀取?32?bits?的數(shù)據(jù)，那么第一個(gè)?thread?讀取的地址，必須是?16*4?=?64?bytes?的倍數(shù)。
如果有一部份的?thread?沒(méi)有讀取內(nèi)存，并不會(huì)影響到其它的?thread?速行?coalesced?的存取。例如：????if(tid?!=?3)?{
????????int?number?=?data[tid];
????}雖然?thread?3?并沒(méi)有讀取數(shù)據(jù)，但是由于其它的?thread?仍符合?coalesced?的條件（假設(shè)?data?的地址是?64?bytes?的倍數(shù)），這樣的內(nèi)存讀取仍會(huì)符合?coalesced?的條件。在目前的?CUDA?1.1?裝置中，每個(gè)?thread?一次讀取的內(nèi)存數(shù)據(jù)量，可以是?32?bits、64?bits、或?128?bits。不過(guò)，32?bits?的效率是最好的。64?bits?的效率會(huì)稍差，而一次讀取?128?bits?的效率則比一次讀取?32?bits?要顯著來(lái)得低（但仍比?non-coalesced?的存取要好）。如果每個(gè)?thread?一次存取的數(shù)據(jù)并不是?32?bits、64?bits、或?128?bits，那就無(wú)法符合?coalesced?的條件。例如，以下的程序：????struct?vec3d?{?float?x,?y,?z;?};
????...
????__global__?void?func(struct?vec3d*?data,?float*?output)
????{
????????output[tid]?=?data[tid].x?*?data[tid].x?+
????????????data[tid].y?*?data[tid].y?+
????????????data[tid].z?*?data[tid].z;
????}
并不是?coalesced?的讀取，因?yàn)?vec3d?的大小是?12?bytes，而非?4?bytes、8?bytes、或?16?bytes。要解決這個(gè)問(wèn)題，可以使用?__align(n)__?的指示，例如：????struct?__align__(16)?vec3d?{?float?x,?y,?z;?};這會(huì)讓?compiler?在?vec3d?后面加上一個(gè)空的?4?bytes，以補(bǔ)齊?16?bytes。另一個(gè)方法，是把數(shù)據(jù)結(jié)構(gòu)轉(zhuǎn)換成三個(gè)連續(xù)的數(shù)組，例如：????__global__?void?func(float*?x,?float*?y,?float*?z,?float*?output)
????{
????????output[tid]?=?x[tid]?*?x[tid]?+?y[tid]?*?y[tid]?+
????????????z[tid]?*?z[tid];
????}如果因?yàn)槠渌蚴箶?shù)據(jù)結(jié)構(gòu)無(wú)法這樣調(diào)整，也可以考慮利用?shared?memory?在?GPU?上做結(jié)構(gòu)的調(diào)整。例如：????__global__?void?func(struct?vec3d*?data,?float*?output)
????{
????????__shared__?float?temp[THREAD_NUM?*?3];
????????const?float*?fdata?=?(float*)?data;
????????temp[tid]?=?fdata[tid];
????????temp[tid?+?THREAD_NUM]?=?fdata[tid?+?THREAD_NUM];
????????temp[tid?+?THREAD_NUM*2]?=?fdata[tid?+?THREAD_NUM*2];
????????__syncthreads();
????????output[tid]?=?temp[tid*3]?*?temp[tid*3]?+
????????????temp[tid*3+1]?*?temp[tid*3+1]?+
????????????temp[tid*3+2]?*?temp[tid*3+2];
????}在上面的例子中，我們先用連續(xù)的方式，把數(shù)據(jù)從?global?memory?讀到?shared?memory。由于?shared?memory?不需要擔(dān)心存取順序（但要注意?bank?conflict?問(wèn)題，參照前一節(jié)），所以可以避開(kāi)?non-coalesced?讀取的問(wèn)題。
CUDA?支援?texture。在?CUDA?的?kernel?程序中，可以利用顯示芯片的?texture?單元，讀取?texture?的數(shù)據(jù)。使用?texture?和?global?memory?最大的差別在于?texture?只能讀取，不能寫入，而且顯示芯片上有一定大小的?texture?cache。因此，讀取?texture?的時(shí)候，不需要符合?coalesced?的規(guī)則，也可以達(dá)到不錯(cuò)的效率。此外，讀取?texture?時(shí)，也可以利用顯示芯片中的?texture?filtering?功能（例如?bilinear?filtering），也可以快速轉(zhuǎn)換數(shù)據(jù)型態(tài)，例如可以直接將?32?bits?RGBA?的數(shù)據(jù)轉(zhuǎn)換成四個(gè)?32?bits?浮點(diǎn)數(shù)。
顯示芯片上的?texture?cache?是針對(duì)一般繪圖應(yīng)用所設(shè)計(jì)，因此它仍最適合有區(qū)塊性質(zhì)的存取動(dòng)作，而非隨機(jī)的存取。因此，同一個(gè)?warp?中的各個(gè)?thread?最好是讀取地址相近的數(shù)據(jù)，才能達(dá)到最高的效率。對(duì)于已經(jīng)能符合?coalesced?規(guī)則的數(shù)據(jù)，使用?global?memory?通常會(huì)比使用?texture?要來(lái)得快。
Stream?processor?里的運(yùn)算單元，基本上是一個(gè)浮點(diǎn)數(shù)的?fused?multiply-add?單元，也就是說(shuō)它可以進(jìn)行一次乘法和一次加法，如下所示：????a?=?b?*?c?+?d;compiler?會(huì)自動(dòng)把適當(dāng)?shù)募臃ê统朔ㄟ\(yùn)算，結(jié)合成一個(gè)?fmad?指令。除了浮點(diǎn)數(shù)的加法及乘法之外，整數(shù)的加法、位運(yùn)算、比較、取最小值、取最大值、及以型態(tài)的轉(zhuǎn)換（浮點(diǎn)數(shù)轉(zhuǎn)整數(shù)或整數(shù)轉(zhuǎn)浮點(diǎn)數(shù)）都是可以全速進(jìn)行的。?整數(shù)的乘法則無(wú)法全速進(jìn)行，但?24?bits?的乘法則可以。在?CUDA?中可以利用內(nèi)建的?__mul24?和?__umul24?函式來(lái)進(jìn)行?24?bits?的整數(shù)乘法。浮點(diǎn)數(shù)的除法是利用先取倒數(shù)，再相乘的方式計(jì)算，因此精確度并不能達(dá)到?IEEE?754?的規(guī)范（最大誤差為?2?ulp）。內(nèi)建的?__fdividef(x,y)?提供更快速的除法，和一般的除法有相同的精確度，但是在?2216?<?y?<?2218?時(shí)會(huì)得到錯(cuò)誤的結(jié)果。此外?CUDA?還提供了一些精確度較低的內(nèi)部函數(shù)，包括?__expf、__logf、__sinf、__cosf、__powf?等等。這些函式的速度較快，但精確度不如標(biāo)準(zhǔn)的函式。詳細(xì)的數(shù)據(jù)可以參考?CUDA?Programming?Guide?1.1?的?Appendix?B。
在?CUDA?中，GPU?不能直接存取主內(nèi)存，只能存取顯卡上的顯示內(nèi)存。因此，會(huì)需要將數(shù)據(jù)從主內(nèi)存先復(fù)制到顯卡內(nèi)存中，進(jìn)行運(yùn)算后，再將結(jié)果從顯卡內(nèi)存中復(fù)制到主內(nèi)存中。這些?復(fù)制的動(dòng)作會(huì)限于?PCI?Express?的速度。使用?PCI?Express?x16?時(shí)，PCI?Express?1.0?可以提供雙向各?4GB/s?的帶寬，而?PCI?Express?2.0?則可提供?8GB/s?的帶寬。當(dāng)然這都是理論值。從一般的內(nèi)存復(fù)制數(shù)據(jù)到顯卡內(nèi)存的時(shí)候，由于一般的內(nèi)存可能隨時(shí)會(huì)被操作系統(tǒng)搬動(dòng)，因此?CUDA?會(huì)先將數(shù)據(jù)復(fù)制到一塊內(nèi)部的內(nèi)存中，才能利用?DMA?將數(shù)據(jù)復(fù)制到顯卡內(nèi)存中。如果想要避免這個(gè)重復(fù)的復(fù)制動(dòng)作，可以使用?cudaMallocHost?函式，在主內(nèi)存中取得一塊?page?locked?的內(nèi)存。不過(guò)，如果要求太大量的?page?locked?的內(nèi)存，將會(huì)影響到操作系統(tǒng)對(duì)內(nèi)存的管理，可能會(huì)減低系統(tǒng)的效率。
雖然矩陣乘法有點(diǎn)老套，不過(guò)因?yàn)樗喈?dāng)簡(jiǎn)單，而且也可以用來(lái)介紹一些有關(guān)?CUDA?的有趣性質(zhì)。
| 矩陣乘法 |
????????for(j?=?0;?j?<?n;?j++)?{
????????????C[i][j]?=?0;
????????????for(k?=?0;?k?<?n;?k++)?{
????????????????C[i][j]?+=?A[i][k]?*?B[k][j];
????????????}
????????}
????}
一開(kāi)始，我們先準(zhǔn)備好產(chǎn)生數(shù)據(jù)、設(shè)定?CUDA?等等的工作：????int?main()
????{
????????float?*a,?*b,?*c,?*d;
????????int?n?=?1000;
????????if(!InitCUDA())?return?0;
????????a?=?(float*)?malloc(sizeof(float)?*?n?*?n);
????????b?=?(float*)?malloc(sizeof(float)?*?n?*?n);
????????c?=?(float*)?malloc(sizeof(float)?*?n?*?n);
????????d?=?(float*)?malloc(sizeof(float)?*?n?*?n);
????????srand(0);
????????matgen(a,?n,?n);
????????matgen(b,?n,?n);
????????clock_t?time?=?matmultCUDA(a,?n,?b,?n,?c,?n,?n);
????????matmult(a,?n,?b,?n,?d,?n,?n);
????????compare_mat(c,?n,?d,?n,?n);
????????double?sec?=?(double)?time?/?CLOCKS_PER_SEC;
????????printf("Time?used:?%.2f?(%.2lf?GFLOPS)\n",?sec,
????????????2.0?*?n?*?n?*?n?/?(sec?*?1E9));
????????return?0;
????}InitCUDA?函式和第一個(gè)?CUDA?程序一樣，可以直接參考前面的文章。以下是上面用到的一些其它的函式：產(chǎn)生矩陣：????void?matgen(float*?a,?int?lda,?int?n)
????{
????????int?i,?j;
????????for(i?=?0;?i?<?n;?i++)?{
????????????for(j?=?0;?j?<?n;?j++)?{
????????????????a[i?*?lda?+?j]?=?(float)?rand()?/?RAND_MAX?+
?????????????????????(float)?rand()?/?(RAND_MAX?*?RAND_MAX);
????????????}
????????}
????}這個(gè)函式只是利用隨機(jī)數(shù)生成器把矩陣填滿?0?~?1?之間的數(shù)字。特別注意到因?yàn)?C?語(yǔ)言中無(wú)法聲明變動(dòng)大小的二維矩陣，所以我們使用?i?*?lda?+?j?的方式。進(jìn)行矩陣乘法：????void?matmult(const?float*?a,?int?lda,?const?float*?b,?int?ldb,
????????float*?c,?int?ldc,?int?n)
????{
????????int?i,?j,?k;
????????for(i?=?0;?i?<?n;?i++)?{
????????????for(j?=?0;?j?<?n;?j++)?{
????????????????double?t?=?0;
????????????????for(k?=?0;?k?<?n;?k++)?{
????????????????????t?+=?a[i?*?lda?+?k]?*?b[k?*?ldb?+?j];
????????????????}
????????????????c[i?*?ldc?+?j]?=?t;
????????????}
????????}
????}這是以?CPU?進(jìn)行矩陣乘法、用來(lái)進(jìn)行驗(yàn)證答案正確與否的程序。特別注意到它用?double?來(lái)儲(chǔ)存暫時(shí)的計(jì)算結(jié)果，以提高精確度。
驗(yàn)證結(jié)果：????void?compare_mat(const?float*?a,?int?lda,
????????const?float*?b,?int?ldb,?int?n)
????{
????????float?max_err?=?0;
????????float?average_err?=?0;
????????int?i,?j;
????????for(i?=?0;?i?<?n;?i++)?{
????????????for(j?=?0;?j?<?n;?j++)?{
????????????????if(b[i?*?ldb?+?j]?!=?0)?{
????????????????????float?err?=?fabs((a[i?*?lda?+?j]?-
?????????????????????????b[i?*?ldb?+?j])?/?b[i?*?ldb?+?j]);
????????????????????if(max_err?<?err)?max_err?=?err;
????????????????????average_err?+=?err;
????????????????}
????????????}
????????}
????????printf("Max?error:?%g?Average?error:?%g\n",
????????????max_err,?average_err?/?(n?*?n));
????}這個(gè)函式計(jì)算兩個(gè)矩陣的最大相對(duì)誤差和平均相對(duì)誤差，并把結(jié)果印出來(lái)。最后是?CUDA?的矩陣乘法的部份：????#define?NUM_THREADS?256
????clock_t?matmultCUDA(const?float*?a,?int?lda,
????????const?float*?b,?int?ldb,?float*?c,?int?ldc,?int?n)
????{
????????float?*ac,?*bc,?*cc;
????????clock_t?start,?end;
????????start?=?clock();
????????cudaMalloc((void**)?&ac,?sizeof(float)?*?n?*?n);
?????????cudaMalloc((void**)?&bc,?sizeof(float)?*?n?*?n);
????????cudaMalloc((void**)?&cc,?sizeof(float)?*?n?*?n);
????????cudaMemcpy2D(ac,?sizeof(float)?*?n,?a,?sizeof(float)?*?lda,
????????????sizeof(float)?*?n,?n,?cudaMemcpyHostToDevice);
????????cudaMemcpy2D(bc,?sizeof(float)?*?n,?b,?sizeof(float)?*?ldb,
????????????sizeof(float)?*?n,?n,?cudaMemcpyHostToDevice);
????????int?blocks?=?(n?+?NUM_THREADS?-?1)?/?NUM_THREADS;
????????matMultCUDA<<<blocks?*?n,?NUM_THREADS>>>
????????????(ac,?n,?bc,?n,?cc,?n,?n);
????????cudaMemcpy2D(c,?sizeof(float)?*?ldc,?cc,?sizeof(float)?*?n,
????????sizeof(float)?*?n,?n,?cudaMemcpyDeviceToHost);
????????cudaFree(ac);
????????cudaFree(bc);
????????cudaFree(cc);
????????end?=?clock();
????????return?end?-?start;
????}這個(gè)函式相當(dāng)單純，就是在顯卡內(nèi)存中配置存放矩陣的內(nèi)存，然后把主內(nèi)存中的矩陣數(shù)據(jù)復(fù)制到顯卡內(nèi)存上。不過(guò)，因?yàn)槲覀兊木仃嚦朔ê娇梢灾付?pitch（即?lda、ldb、和?ldc），所以如果用一般的?cudaMemcpy?函式來(lái)復(fù)制內(nèi)存的話，會(huì)需要每個(gè)?row?都分開(kāi)復(fù)制，那會(huì)需要呼叫很多次?cudaMemcpy?函式，會(huì)使效率變得很差。因此，在這里我們用了一個(gè)新的?cudaMemcpy2D?函式，它是用來(lái)復(fù)制二維數(shù)組，可以指定數(shù)組的?pitch。這樣就可以透過(guò)一次函數(shù)調(diào)用就可以了。進(jìn)行計(jì)算的?kernel?如下：????__global__?static?void?matMultCUDA(const?float*?a,?size_t?lda,
????????const?float*?b,?size_t?ldb,?float*?c,?size_t?ldc,?int?n)
????{
????????const?int?tid?=?threadIdx.x;
????????const?int?bid?=?blockIdx.x;
????????const?int?idx?=?bid?*?blockDim.x?+?tid;
????????const?int?row?=?idx?/?n;
????????const?int?column?=?idx?%?n;
????????int?i;
????????if(row?<?n?&&?column?<?n)?{
????????????float?t?=?0;
????????????for(i?=?0;?i?<?n;?i++)?{
????????????????t?+=?a[row?*?lda?+?i]?*?b[i?*?ldb?+?column];
????????????}
????????????c[row?*?ldc?+?column]?=?t;
????????}
????}這個(gè)函式一開(kāi)始先從?bid?和?tid?計(jì)算出這個(gè)?thread?應(yīng)該計(jì)算的?row?和?column，在判斷?row?和?column?在范圍內(nèi)之后，就直接進(jìn)行計(jì)算，并把結(jié)果寫到?c?矩陣中，是非常單純的函式。在?GeForce?8800GT?上實(shí)際執(zhí)行的結(jié)果如下：????Max?error:?2.01484e-006?Average?error:?3.36637e-007
????Time?used:?1.1560?(1.73?GFLOPS)可以看到兩個(gè)問(wèn)題：1?很明顯的，執(zhí)行效率相當(dāng)?shù)吐洹?/span>2?最大相對(duì)誤差偏高。理想上應(yīng)該要低于?1e-6。計(jì)算結(jié)果的誤差偏高的原因是，在?CPU?上進(jìn)行計(jì)算時(shí)，我們使用?double（即?64?bits?浮點(diǎn)數(shù)）來(lái)累進(jìn)計(jì)算過(guò)程，而在?GPU?上則只能用?float（32?bits?浮點(diǎn)數(shù)）。在累加大量數(shù)字的時(shí)候，由于累加結(jié)果很快會(huì)變大，因此后面的數(shù)字很容易被舍去過(guò)多的位數(shù)。由于?CUDA?的浮點(diǎn)數(shù)運(yùn)算，在進(jìn)行加、減、乘法時(shí)是符合?IEEE?754?規(guī)定的精確度的，因此，我們可以利用?Kahan's?Summation?Formula?來(lái)提高精確度。把程序改成：????if(row?<?n?&&?column?<?n)?{
????????float?t?=?0;
????????float?y?=?0;
????????for(i?=?0;?i?<?n;?i++)?{
????????????float?r;
????????????y?-=?a[row?*?lda?+?i]?*?b[i?*?ldb?+?column];
????????????r?=?t?-?y;
????????????y?=?(r?-?t)?+?y;
????????????t?=?r;
????????}
????}修改后的程序的執(zhí)行結(jié)果是：????Max?error:?1.19209e-007?Average?error:?4.22751e-008
????Time?used:?1.1560?(1.73?GFLOPS)可以看到相對(duì)誤差有很大的改善，效率則沒(méi)什么變化。由于?Kahan's?Summation?Formula?需要的運(yùn)算量提高，但是效率卻沒(méi)有什么改變，可以看出這個(gè)?kernel?主要的瓶頸應(yīng)該是在內(nèi)存的存取動(dòng)作上。這是因?yàn)橛写罅康膬?nèi)存讀取是重復(fù)的。例如，矩陣?a?的一個(gè)?row?在每次進(jìn)行計(jì)算時(shí)都被重復(fù)讀入，但這是相當(dāng)浪費(fèi)的。這樣的計(jì)算方式，總共需要讀取?2*n3?次內(nèi)存。如果讓一個(gè)?row?只需要讀入一次的話，就可以減到為?n3+n2?次。
| 第一個(gè)改良 |
????????????(ac,?n,?bc,?n,?cc,?n,?n);kernel?的部份則改成：????__global__?static?void?matMultCUDA(const?float*?a,?size_t?lda,
????????const?float*?b,?size_t?ldb,?float*?c,?size_t?ldc,?int?n)
????{
????????extern?__shared__?float?data[];
????????const?int?tid?=?threadIdx.x;
????????const?int?row?=?blockIdx.x;
????????int?i,?j;
????????for(i?=?tid;?i?<?n;?i?+=?blockDim.x)?{
????????????data[i]?=?a[row?*?lda?+?i];
????????}
????????__syncthreads();
????????for(j?=?tid;?j?<?n;?j?+=?blockDim.x)?{
????????????float?t?=?0;
????????????float?y?=?0;
????????????for(i?=?0;?i?<?n;?i++)?{
????????????????float?r;
????????????????y?-=?data[i]?*?b[i?*?ldb?+?j];
????????????????r?=?t?-?y;
????????????????y?=?(r?-?t)?+?y;
????????????????t?=?r;
????????????}
????????????c[row?*?ldc?+?j]?=?t;
????????}
????}
第一個(gè)部份先把整個(gè)?row?讀到?shared?memory?中，而第二個(gè)部份則進(jìn)行計(jì)算，并沒(méi)有太大的變化。主要的差別是現(xiàn)在一個(gè)?row?只由一個(gè)?block?進(jìn)行計(jì)算。在?GeForce?8800GT?上，執(zhí)行的結(jié)果是：????Max?error:?1.19209e-007??Average?error:?4.22751e-008
????Time?used:?0.4220???(4.74?GFLOPS)很明顯的，計(jì)算的結(jié)果并沒(méi)有改變，不過(guò)速度則提高了超過(guò)一倍。雖然如此，但是這樣的效率仍不盡理想，因?yàn)槔碚撋?GeForce?8800GT?有超過(guò)?300GFLOPS?的運(yùn)算性能。即使是把?Kahan's?Summation?Formula?所需要的額外運(yùn)算考慮進(jìn)去，這樣的效率仍然連理論最大值的十分之一都不到。會(huì)有這樣的結(jié)果，原因其實(shí)還是同樣的：對(duì)內(nèi)存的存取次數(shù)太多了。雖然現(xiàn)在?A?矩陣的?row?的數(shù)據(jù)已經(jīng)不再需要重復(fù)讀取，但是?B?矩陣的?column?的數(shù)據(jù)仍然一直被重復(fù)讀取。另一個(gè)問(wèn)題比較不是那么明顯：對(duì)?B?矩陣的讀取，雖然看起來(lái)不連續(xù)，但實(shí)際上它是連續(xù)的。這是因?yàn)椴煌?thread?會(huì)讀取不同的?column，因此同時(shí)間每個(gè)?thread?讀取的各個(gè)?column?加起來(lái)，就是一個(gè)連續(xù)的內(nèi)存區(qū)塊。那么，為什么效率還是不佳呢？這是因?yàn)?#xff0c;GPU?上的內(nèi)存控制器，從某個(gè)固定的倍數(shù)地址開(kāi)始讀取，才會(huì)有最高的效率（例如?16?bytes?的倍數(shù)）。由于矩陣大小并不是?16?的倍數(shù)（這里使用的是?1000x1000?的矩陣），所以造成效率不佳的情形。要解決這個(gè)問(wèn)題，我們可以在?cudaMalloc的時(shí)候稍微修改一下，讓寬度變成?適當(dāng)?shù)谋稊?shù)就可以了。但是，適當(dāng)?shù)谋稊?shù)是多少呢？幸運(yùn)的是，我們并不需要知道這些細(xì)節(jié)。CUDA?提供了一個(gè)?cudaMallocPitch的函式，可以自動(dòng)以最佳的倍數(shù)來(lái)配置內(nèi)存。因此，我們可以把?cudaMalloc?的部份改成：????size_t?pitch_a,?pitch_b,?pitch_c;
????cudaMallocPitch((void**)?&ac,?&pitch_a,?sizeof(float)?*?n,?n);
????cudaMallocPitch((void**)?&bc,?&pitch_b,?sizeof(float)?*?n,?n);
????cudaMallocPitch((void**)?&cc,?&pitch_c,?sizeof(float)?*?n,?n);cudaMallocPitch?函式會(huì)以適當(dāng)?shù)谋稊?shù)配置內(nèi)存，并把配置的寬度傳回。因此，在把矩陣復(fù)制到顯卡內(nèi)存上時(shí)，要使用它傳回的寬度：????cudaMemcpy2D(ac,?pitch_a,?a,?sizeof(float)?*?lda,
????????sizeof(float)?*?n,?n,?cudaMemcpyHostToDevice);
????cudaMemcpy2D(bc,?pitch_b,?b,?sizeof(float)?*?ldb,
????????sizeof(float)?*?n,?n,?cudaMemcpyHostToDevice);呼叫?kernel?的部份也需要修改：????matMultCUDA<<<n,?NUM_THREADS,?sizeof(float)?*?n>>>
????????(ac,?pitch_a?/?sizeof(float),?bc,?pitch_b?/?sizeof(float),
????????cc,?pitch_c?/?sizeof(float),?n);同樣的，把計(jì)算結(jié)果復(fù)制回到主內(nèi)存時(shí)，也要使用傳回的寬度值：????cudaMemcpy2D(c,?sizeof(float)?*?ldc,?cc,?pitch_c,
????????sizeof(float)?*?n,?n,?cudaMemcpyDeviceToHost);這樣就完成了。Kernel?部份則不需要修改。這樣的修改有多大的效果呢？在?GeForce?8800GT?上執(zhí)行，結(jié)果如下：????Max?error:?1.19209e-007??Average?error:?4.22751e-008
????Time?used:?0.1250???(16.00?GFLOPS)可以看到，執(zhí)行速度又再大幅提高了三倍多！而這只是把內(nèi)存的配置方式稍微修改一下而已。雖然執(zhí)行速度提高了很多，但是，和前面提到的理論值相比，其實(shí)還是有相當(dāng)?shù)牟罹?。這是因?yàn)?#xff0c;前面也提到過(guò)，這樣的做法需要?n3+n2次的內(nèi)存讀取，和?n2次?的內(nèi)存寫入動(dòng)作。由于?n?=?1000，每個(gè)數(shù)字的大小是?32?bits，所以總共的內(nèi)存存取數(shù)據(jù)量約為?4GB。除以實(shí)際執(zhí)行的時(shí)間?0.125?秒，得到的帶寬數(shù)值是約?32GB/s，這已經(jīng)接近?GeForce?8800GT?顯卡內(nèi)存的帶寬了。由于我們計(jì)算時(shí)間的時(shí)候，把配置內(nèi)存、以及數(shù)據(jù)的復(fù)制動(dòng)作也計(jì)算進(jìn)去，因此實(shí)際上花費(fèi)在?kernel?的時(shí)間是更短的（約?0.09?秒）。因此，可以很明顯的看出，這個(gè)程序的效率，是受限于內(nèi)存帶寬的。
| 進(jìn)一步的改良 |
????+?A(0~15,?32~47)?*?B(32~47,?0~15)?+?...這樣一來(lái)，我們就可以把兩個(gè)小矩陣加載到?shared?memory，則小矩陣本身的乘法就不需要再存取任何外部的內(nèi)存了！這樣一來(lái)，假設(shè)小矩陣的大小是?k，則實(shí)際上需要的內(nèi)存存取次數(shù)就會(huì)變成約?2k2(n/k)3?=?2n3/k。由于目前?CUDA?每個(gè)?block?的?thread?數(shù)目最多是?512，因此?k?=?16?似乎是一個(gè)相當(dāng)理想的數(shù)字（共?256?個(gè)?threads）。因此，對(duì)于一個(gè)?n?=?1000?的矩陣來(lái)說(shuō)，我們可以把內(nèi)存存取的量減少到約?500MB，也就是上一節(jié)的存取量的?1/8。理論上，這樣應(yīng)該可以讓效率提高八倍（假設(shè)沒(méi)有遇到別的瓶頸）。為了方便進(jìn)行區(qū)塊的計(jì)算，我們讓每個(gè)?block?有?16x16?個(gè)?threads，再建立?(n/16)x(n/16)?個(gè)?blocks。把呼叫?kernel?的地方改成：????int?bx?=?(n?+?BLOCK_SIZE?-?1)?/?BLOCK_SIZE;
????dim3?blocks(bx,?bx);
????dim3?threads(BLOCK_SIZE,?BLOCK_SIZE);
????matMultCUDA<<<blocks,?threads>>>(ac,?pitch_a?/?sizeof(float),
????????bc,?pitch_b?/?sizeof(float),?cc,?pitch_c?/?sizeof(float),?n);BLOCK_SIZE?則是定義成?16。dim3?是?CUDA?的一種數(shù)據(jù)型態(tài)，表示一個(gè)?3D?的向量。在這里，我們透過(guò)?dim3?來(lái)建立?16x16?個(gè)?threads?的?block，和?(n/16)x(n/16)?個(gè)?blocks。Kernel?程序的部份，則改成：????__global__?static?void?matMultCUDA(const?float*?a,?size_t?lda,
????????const?float*?b,?size_t?ldb,?float*?c,?size_t?ldc,?int?n)
????{
????????__shared__?float?matA[BLOCK_SIZE][BLOCK_SIZE];
????????__shared__?float?matB[BLOCK_SIZE][BLOCK_SIZE];
????????const?int?tidc?=?threadIdx.x;
????????const?int?tidr?=?threadIdx.y;
????????const?int?bidc?=?blockIdx.x?*?BLOCK_SIZE;
????????const?int?bidr?=?blockIdx.y?*?BLOCK_SIZE;
????????int?i,?j;
????????float?results?=?0;
????????float?comp?=?0;
????????for(j?=?0;?j?<?n;?j?+=?BLOCK_SIZE)?{
????????????if(tidr?+?bidr?<?n?&&?tidc?+?j?<?n)?{
????????????????matA[tidr][tidc]?=?a[(tidr?+?bidr)?*?lda?+?tidc?+?j];
????????????}
????????????else?{
????????????????matA[tidr][tidc]?=?0;
????????????}
????????????if(tidr?+?j?<?n?&&?tidc?+?bidc?<?n)?{
????????????????matB[tidr][tidc]?=?b[(tidr?+?j)?*?ldb?+?tidc?+?bidc];
????????????}
????????????else?{
????????????????matB[tidr][tidc]?=?0;
????????????}
????????????__syncthreads();
????????????for(i?=?0;?i?<?BLOCK_SIZE;?i++)?{
????????????????float?t;
????????????????comp?-=?matA[tidr][i]?*?matB[i][tidc];
????????????????t?=?results?-?comp;
????????????????comp?=?(t?-?results)?+?comp;
????????????????results?=?t;
????????????}
????????????__syncthreads();
????????}
????????if(tidr?+?bidr?<?n?&&?tidc?+?bidc?<?n)?{
????????????c[(tidr?+?bidr)?*?ldc?+?tidc?+?bidc]?=?results;
????????}
????}注意到因?yàn)槲覀儸F(xiàn)在使用?16x16?的?threads，因此?threadIdx?變量可以取得?threadIdx.x?和?threadIdx.y，范圍分別是?0?~?15。blockIdx.x?和?blockIdx.y?變量也是同樣的情形，范圍分別是?0?~?n/16。在程序中，因?yàn)榫仃嚨拇笮〔灰欢〞?huì)是?16?的倍數(shù)，因此需要使用?if?判斷式檢查是否超出矩陣范圍。這個(gè)版本在?GeForce?8800GT?上的執(zhí)行結(jié)果如下：????Max?error:?1.19209e-007??Average?error:?4.22751e-008
????Time?used:?0.0780???(25.64?GFLOPS)速度雖然提高了，但是似乎并沒(méi)有達(dá)到預(yù)期中的八倍。當(dāng)然，前面提到過(guò)，我們?cè)谟?jì)算時(shí)間時(shí)，把一些復(fù)制內(nèi)存、配置內(nèi)存的動(dòng)作也計(jì)算在內(nèi)，這些動(dòng)作的時(shí)?間并不會(huì)縮短。實(shí)際上?kernel?的運(yùn)行時(shí)間，大約是?0.053?秒左右（約略相當(dāng)于?38GFLOPS），比上一節(jié)的版本快了將近一倍。如果這一版程序已經(jīng)不再限于內(nèi)存帶寬，那為什么沒(méi)有達(dá)到預(yù)期的效率呢？這是因?yàn)檫@一版程序已經(jīng)是限于指令周期了。除了使用?Kahan's?Summation?Formula?會(huì)需要更多的運(yùn)算之外，程序中也有大量計(jì)算矩陣地址的乘法等等，這都會(huì)需要花費(fèi)運(yùn)算資源。另外，那些用來(lái)判斷超出矩陣范圍的?if?判斷式，也會(huì)有一定的影響。要把那些?if?判斷式去掉，有一個(gè)方法是，在配置內(nèi)存時(shí)，就配置成?16?的倍數(shù)，并在復(fù)制矩陣到顯卡內(nèi)存之前，先將它清為?0。如下所示：????int?newn?=?((n?+?BLOCK_SIZE?-?1)?/?BLOCK_SIZE)?*?BLOCK_SIZE;
????cudaMallocPitch((void**)?&ac,?&pitch_a,
?????????sizeof(float)?*?newn,?newn);
????cudaMallocPitch((void**)?&bc,?&pitch_b,
????????sizeof(float)?*?newn,?newn);
????cudaMallocPitch((void**)?&cc,?&pitch_c,
????????sizeof(float)?*?newn,?newn);
????cudaMemset(ac,?0,?pitch_a?*?newn);
????cudaMemset(bc,?0,?pitch_b?*?newn);?這樣一來(lái)，我們就可以把?kernel?中的?if?判斷式都移除了：????__global__?static?void?matMultCUDA(const?float*?a,?size_t?lda,
????????const?float*?b,?size_t?ldb,?float*?c,?size_t?ldc,?int?n)
????{
????????__shared__?float?matA[BLOCK_SIZE][BLOCK_SIZE];
????????__shared__?float?matB[BLOCK_SIZE][BLOCK_SIZE];
????????const?int?tidc?=?threadIdx.x;
????????const?int?tidr?=?threadIdx.y;
????????const?int?bidc?=?blockIdx.x?*?BLOCK_SIZE;
????????const?int?bidr?=?blockIdx.y?*?BLOCK_SIZE;
????????int?i,?j;
????????float?results?=?0;
????????float?comp?=?0;
????????for(j?=?0;?j?<?n;?j?+=?BLOCK_SIZE)?{
????????????matA[tidr][tidc]?=?a[(tidr?+?bidr)?*?lda?+?tidc?+?j];
????????????matB[tidr][tidc]?=?b[(tidr?+?j)?*?ldb?+?tidc?+?bidc];
????????????__syncthreads();
????????????for(i?=?0;?i?<?BLOCK_SIZE;?i++)?{
????????????????float?t;
????????????????comp?-=?matA[tidr][i]?*?matB[i][tidc];
????????????????t?=?results?-?comp;
????????????????comp?=?(t?-?results)?+?comp;
????????????????results?=?t;
????????????}
????????????__syncthreads();
????????}
????????c[(tidr?+?bidr)?*?ldc?+?tidc?+?bidc]?=?results;
????}這個(gè)版本的執(zhí)行結(jié)果是：????Max?error:?1.19209e-007??Average?error:?4.22751e-008
????Time?used:?0.0780???(25.64?GFLOPS)似乎沒(méi)有改善。不過(guò)，實(shí)際上?kernel?的運(yùn)行時(shí)間已經(jīng)減少到?0.042?秒（約略相當(dāng)于?48GFLOPS）。
| 結(jié)論 |
| GPU?的基本介紹 |
詳細(xì)的?multiprocessor?信息，都可以透過(guò)?CUDA?的?cudaGetDeviceProperties()?函式或?cuDeviceGetProperties()?函式取得。不過(guò)，目前還沒(méi)有辦法直接取得一個(gè)顯示芯片中有多少?multiprocessor?的信息。
在?CUDA?中，大部份基本的運(yùn)算動(dòng)作，都可以由?stream?processor?進(jìn)行。每個(gè)?stream?processor?都包含一個(gè)?FMA（fused-multiply-add）單元，可以進(jìn)行一個(gè)乘法和一個(gè)加法。比較復(fù)雜的運(yùn)算則會(huì)需要比較長(zhǎng)的時(shí)間。
| 執(zhí)行過(guò)程 |
實(shí)際上，雖然一個(gè)?multiprocessor?只有八個(gè)?stream?processor，但是由于?stream?processor?進(jìn)行各種運(yùn)算都有?latency，更不用提內(nèi)存存取的?latency，因此?CUDA?在執(zhí)行程序的時(shí)候，是以?warp?為單位。目前的?CUDA?裝置，一個(gè)?warp?里面有?32?個(gè)?threads，分成兩組?16?threads?的?half-warp。由于?stream?processor?的運(yùn)算至少有?4?cycles?的?latency，因此對(duì)一個(gè)?4D?的?stream?processors?來(lái)說(shuō)，一次至少執(zhí)行?16?個(gè)?threads（即?half-warp）才能有效隱藏各種運(yùn)算的?latency。
由于?multiprocessor?中并沒(méi)有太多別的內(nèi)存，因此每個(gè)?thread?的狀態(tài)都是直接保存在?multiprocessor?的寄存器中。所以，如果一個(gè)?multiprocessor?同時(shí)有愈多的?thread?要執(zhí)行，就會(huì)需要愈多的寄存器空間。例如，假設(shè)一個(gè)?block?里面有?256?個(gè)?threads，每個(gè)?thread?用到?20?個(gè)寄存器，那么總共就需要?256x20?=?5,120?個(gè)寄存器才能保存每個(gè)?thread?的狀態(tài)。目前?CUDA?裝置中每個(gè)?multiprocessor?有?8,192?個(gè)寄存器，因此，如果每個(gè)?thread?使用到?16?個(gè)寄存器，那就表示一個(gè)?multiprocessor?同時(shí)最多只能維持?512?個(gè)?thread?的執(zhí)行。如果同時(shí)進(jìn)行的?thread?數(shù)目超過(guò)這個(gè)數(shù)字，那么就會(huì)需要把一部份的數(shù)據(jù)儲(chǔ)存在顯卡內(nèi)存中，就會(huì)降低執(zhí)行的效率了。編者注：在NVIDIA?GT200中的Register?File大小增加了一倍，在FP32下可用的register?file為16K，FP64下是8K。
| Shared?memory |
????...
????int?number?=?data[16?*?tid];
會(huì)造成嚴(yán)重的?bank?conflict，為了避免這個(gè)問(wèn)題，可以把數(shù)據(jù)的排列方式稍加修改，把存取方式改成：????int?row?=?tid?/?16;
????int?column?=?tid?%?16;
????data[row?*?17?+?column]?=?global_data[tid];
????...
????int?number?=?data[17?*?tid];這樣就不會(huì)造成?bank?conflict?了。編者注：share?memory在NVIDIA的文檔中其實(shí)還有不同的叫法，例如PDC（Parallel?Data?Cache）、PBSM（per-block?share?memory）。
| Global?memory |
如果有一部份的?thread?沒(méi)有讀取內(nèi)存，并不會(huì)影響到其它的?thread?速行?coalesced?的存取。例如：????if(tid?!=?3)?{
????????int?number?=?data[tid];
????}雖然?thread?3?并沒(méi)有讀取數(shù)據(jù)，但是由于其它的?thread?仍符合?coalesced?的條件（假設(shè)?data?的地址是?64?bytes?的倍數(shù)），這樣的內(nèi)存讀取仍會(huì)符合?coalesced?的條件。在目前的?CUDA?1.1?裝置中，每個(gè)?thread?一次讀取的內(nèi)存數(shù)據(jù)量，可以是?32?bits、64?bits、或?128?bits。不過(guò)，32?bits?的效率是最好的。64?bits?的效率會(huì)稍差，而一次讀取?128?bits?的效率則比一次讀取?32?bits?要顯著來(lái)得低（但仍比?non-coalesced?的存取要好）。如果每個(gè)?thread?一次存取的數(shù)據(jù)并不是?32?bits、64?bits、或?128?bits，那就無(wú)法符合?coalesced?的條件。例如，以下的程序：????struct?vec3d?{?float?x,?y,?z;?};
????...
????__global__?void?func(struct?vec3d*?data,?float*?output)
????{
????????output[tid]?=?data[tid].x?*?data[tid].x?+
????????????data[tid].y?*?data[tid].y?+
????????????data[tid].z?*?data[tid].z;
????}
并不是?coalesced?的讀取，因?yàn)?vec3d?的大小是?12?bytes，而非?4?bytes、8?bytes、或?16?bytes。要解決這個(gè)問(wèn)題，可以使用?__align(n)__?的指示，例如：????struct?__align__(16)?vec3d?{?float?x,?y,?z;?};這會(huì)讓?compiler?在?vec3d?后面加上一個(gè)空的?4?bytes，以補(bǔ)齊?16?bytes。另一個(gè)方法，是把數(shù)據(jù)結(jié)構(gòu)轉(zhuǎn)換成三個(gè)連續(xù)的數(shù)組，例如：????__global__?void?func(float*?x,?float*?y,?float*?z,?float*?output)
????{
????????output[tid]?=?x[tid]?*?x[tid]?+?y[tid]?*?y[tid]?+
????????????z[tid]?*?z[tid];
????}如果因?yàn)槠渌蚴箶?shù)據(jù)結(jié)構(gòu)無(wú)法這樣調(diào)整，也可以考慮利用?shared?memory?在?GPU?上做結(jié)構(gòu)的調(diào)整。例如：????__global__?void?func(struct?vec3d*?data,?float*?output)
????{
????????__shared__?float?temp[THREAD_NUM?*?3];
????????const?float*?fdata?=?(float*)?data;
????????temp[tid]?=?fdata[tid];
????????temp[tid?+?THREAD_NUM]?=?fdata[tid?+?THREAD_NUM];
????????temp[tid?+?THREAD_NUM*2]?=?fdata[tid?+?THREAD_NUM*2];
????????__syncthreads();
????????output[tid]?=?temp[tid*3]?*?temp[tid*3]?+
????????????temp[tid*3+1]?*?temp[tid*3+1]?+
????????????temp[tid*3+2]?*?temp[tid*3+2];
????}在上面的例子中，我們先用連續(xù)的方式，把數(shù)據(jù)從?global?memory?讀到?shared?memory。由于?shared?memory?不需要擔(dān)心存取順序（但要注意?bank?conflict?問(wèn)題，參照前一節(jié)），所以可以避開(kāi)?non-coalesced?讀取的問(wèn)題。
| Texture |
顯示芯片上的?texture?cache?是針對(duì)一般繪圖應(yīng)用所設(shè)計(jì)，因此它仍最適合有區(qū)塊性質(zhì)的存取動(dòng)作，而非隨機(jī)的存取。因此，同一個(gè)?warp?中的各個(gè)?thread?最好是讀取地址相近的數(shù)據(jù)，才能達(dá)到最高的效率。對(duì)于已經(jīng)能符合?coalesced?規(guī)則的數(shù)據(jù)，使用?global?memory?通常會(huì)比使用?texture?要來(lái)得快。
| 運(yùn)算單元 |
| 和主內(nèi)存間的數(shù)據(jù)傳輸 |
轉(zhuǎn)載于:https://blog.51cto.com/390820/178589
總結(jié)
以上是生活随笔為你收集整理的深入浅出谈CUDA（二）的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。
- 上一篇： Sql server 2005 中的de
- 下一篇： System Center Virtua