前面介紹的計算平方和的程序，似乎沒有什么實用價值。所以我們的第二個?CUDA?程序，要做一個確實有（某些）實用價值的程序，也就是進行矩陣乘法。而且，這次我們會使用浮點數。

雖然矩陣乘法有點老套，不過因為它相當簡單，而且也可以用來介紹一些有關?CUDA?的有趣性質。

矩陣乘法

為了單純起見，我們這里以方形的矩陣為例子。基本上，假設有兩個矩陣?A?和?B，則計算?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];
????????????}
????????}
????}

一開始，我們先準備好產生數據、設定?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?函式和第一個?CUDA?程序一樣，可以直接參考前面的文章。以下是上面用到的一些其它的函式：

產生矩陣：

????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);
????????????}
????????}
????}

這個函式只是利用隨機數生成器把矩陣填滿?0?~?1?之間的數字。特別注意到因為?C?語言中無法聲明變動大小的二維矩陣，所以我們使用?i?*?lda?+?j?的方式。

進行矩陣乘法：

????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?進行矩陣乘法、用來進行驗證答案正確與否的程序。特別注意到它用?double?來儲存暫時的計算結果，以提高精確度。

驗證結果：

????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));
????}

這個函式計算兩個矩陣的最大相對誤差和平均相對誤差，并把結果印出來。

最后是?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;
????}

這個函式相當單純，就是在顯卡內存中配置存放矩陣的內存，然后把主內存中的矩陣數據復制到顯卡內存上。不過，因為我們的矩陣乘法函式可以指定?pitch（即?lda、ldb、和?ldc），所以如果用一般的?cudaMemcpy?函式來復制內存的話，會需要每個?row?都分開復制，那會需要呼叫很多次?cudaMemcpy?函式，會使效率變得很差。因此，在這里我們用了一個新的?cudaMemcpy2D?函式，它是用來復制二維數組，可以指定數組的?pitch。這樣就可以透過一次函數調用就可以了。

進行計算的?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;
????????}
????}

這個函式一開始先從?bid?和?tid?計算出這個?thread?應該計算的?row?和?column，在判斷?row?和?column?在范圍內之后，就直接進行計算，并把結果寫到?c?矩陣中，是非常單純的函式。

在?GeForce?8800GT?上實際執行的結果如下：

????Max?error:?2.01484e-006?Average?error:?3.36637e-007
????Time?used:?1.1560?(1.73?GFLOPS)

可以看到兩個問題：

1?很明顯的，執行效率相當低落。

2?最大相對誤差偏高。理想上應該要低于?1e-6。

計算結果的誤差偏高的原因是，在?CPU?上進行計算時，我們使用?double（即?64?bits?浮點數）來累進計算過程，而在?GPU?上則只能用?float（32?bits?浮點數）。在累加大量數字的時候，由于累加結果很快會變大，因此后面的數字很容易被舍去過多的位數。

由于?CUDA?的浮點數運算，在進行加、減、乘法時是符合?IEEE?754?規定的精確度的，因此，我們可以利用?Kahan's?Summation?Formula?來提高精確度。把程序改成：

????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;
????????}
????}

修改后的程序的執行結果是：

????Max?error:?1.19209e-007?Average?error:?4.22751e-008
????Time?used:?1.1560?(1.73?GFLOPS)

可以看到相對誤差有很大的改善，效率則沒什么變化。

由于?Kahan's?Summation?Formula?需要的運算量提高，但是效率卻沒有什么改變，可以看出這個?kernel?主要的瓶頸應該是在內存的存取動作上。這是因為有大量的內存讀取是重復的。例如，矩陣?a?的一個?row?在每次進行計算時都被重復讀入，但這是相當浪費的。這樣的計算方式，總共需要讀取?2*n3?次內存。如果讓一個?row?只需要讀入一次的話，就可以減到為?n3+n2?次。

第一個改良

?

和我們的第一個?CUDA?程序一樣，我們可以利用?shared?memory?來儲存每個?row?的數據。不過，因為只有同一個?block?的?threads?可以共享?shared?memory，因此現在一個?row?只能由同一個?block?的?threads?來進行計算。另外我們也需要能存放一整個?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;
????????}
????}

第一個部份先把整個?row?讀到?shared?memory?中，而第二個部份則進行計算，并沒有太大的變化。主要的差別是現在一個?row?只由一個?block?進行計算。

在?GeForce?8800GT?上，執行的結果是：

????Max?error:?1.19209e-007??Average?error:?4.22751e-008
????Time?used:?0.4220???(4.74?GFLOPS)

很明顯的，計算的結果并沒有改變，不過速度則提高了超過一倍。雖然如此，但是這樣的效率仍不盡理想，因為理論上?GeForce?8800GT?有超過?300GFLOPS?的運算性能。即使是把?Kahan's?Summation?Formula?所需要的額外運算考慮進去，這樣的效率仍然連理論最大值的十分之一都不到。

會有這樣的結果，原因其實還是同樣的：對內存的存取次數太多了。雖然現在?A?矩陣的?row?的數據已經不再需要重復讀取，但是?B?矩陣的?column?的數據仍然一直被重復讀取。

另一個問題比較不是那么明顯：對?B?矩陣的讀取，雖然看起來不連續，但實際上它是連續的。這是因為不同的?thread?會讀取不同的?column，因此同時間每個?thread?讀取的各個?column?加起來，就是一個連續的內存區塊。那么，為什么效率還是不佳呢？這是因為，GPU?上的內存控制器，從某個固定的倍數地址開始讀取，才會有最高的效率（例如?16?bytes?的倍數）。由于矩陣大小并不是?16?的倍數（這里使用的是?1000x1000?的矩陣），所以造成效率不佳的情形。

要解決這個問題，我們可以在?cudaMalloc的時候稍微修改一下，讓寬度變成?適當的倍數就可以了。但是，適當的倍數是多少呢？幸運的是，我們并不需要知道這些細節。CUDA?提供了一個?cudaMallocPitch的函式，可以自動以最佳的倍數來配置內存。因此，我們可以把?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?函式會以適當的倍數配置內存，并把配置的寬度傳回。因此，在把矩陣復制到顯卡內存上時，要使用它傳回的寬度：

????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);

同樣的，把計算結果復制回到主內存時，也要使用傳回的寬度值：

????cudaMemcpy2D(c,?sizeof(float)?*?ldc,?cc,?pitch_c,
????????sizeof(float)?*?n,?n,?cudaMemcpyDeviceToHost);

這樣就完成了。Kernel?部份則不需要修改。

這樣的修改有多大的效果呢？在?GeForce?8800GT?上執行，結果如下：

????Max?error:?1.19209e-007??Average?error:?4.22751e-008
????Time?used:?0.1250???(16.00?GFLOPS)

可以看到，執行速度又再大幅提高了三倍多！而這只是把內存的配置方式稍微修改一下而已。

雖然執行速度提高了很多，但是，和前面提到的理論值相比，其實還是有相當的差距。這是因為，前面也提到過，這樣的做法需要?n3+n2次的內存讀取，和?n2次?的內存寫入動作。由于?n?=?1000，每個數字的大小是?32?bits，所以總共的內存存取數據量約為?4GB。除以實際執行的時間?0.125?秒，得到的帶寬數值是約?32GB/s，這已經接近?GeForce?8800GT?顯卡內存的帶寬了。由于我們計算時間的時候，把配置內存、以及數據的復制動作也計算進去，因此實際上花費在?kernel?的時間是更短的（約?0.09?秒）。因此，可以很明顯的看出，這個程序的效率，是受限于內存帶寬的。

進一步的改良

上一節的結論顯示出，矩陣乘法的程序，效率是受限于內存帶寬的。那有沒有辦法降低內存的存取次數呢？答案當然是有的，不然就不會有這一節了?:)

要進一步降低內存帶寬的使用，可以注意到，在上一節的方法中，雖然?A?矩陣的存取次數被減至最低，但是?B?矩陣的存取次數并沒有減少。這是因為我們只將?A?矩陣的?row?加載到?shared?memory?中，但是?B?矩陣的?column?也是有被重復使用的。理想上應該也可以避免重復加載才對。不過，由于?B?矩陣的?column?使用的時機，和?A?矩陣的?row?是不同的，所以并不能直接這樣做。

解決方法是?"blocking"。也就是把整個矩陣乘法的動作，切割成很多小矩陣的乘法。例如，要計算?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)?+?...

這樣一來，我們就可以把兩個小矩陣加載到?shared?memory，則小矩陣本身的乘法就不需要再存取任何外部的內存了！這樣一來，假設小矩陣的大小是?k，則實際上需要的內存存取次數就會變成約?2k2(n/k)3?=?2n3/k。

由于目前?CUDA?每個?block?的?thread?數目最多是?512，因此?k?=?16?似乎是一個相當理想的數字（共?256?個?threads）。因此，對于一個?n?=?1000?的矩陣來說，我們可以把內存存取的量減少到約?500MB，也就是上一節的存取量的?1/8。理論上，這樣應該可以讓效率提高八倍（假設沒有遇到別的瓶頸）。

為了方便進行區塊的計算，我們讓每個?block?有?16x16?個?threads，再建立?(n/16)x(n/16)?個?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?的一種數據型態，表示一個?3D?的向量。在這里，我們透過?dim3?來建立?16x16?個?threads?的?block，和?(n/16)x(n/16)?個?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;
????????}
????}

注意到因為我們現在使用?16x16?的?threads，因此?threadIdx?變量可以取得?threadIdx.x?和?threadIdx.y，范圍分別是?0?~?15。blockIdx.x?和?blockIdx.y?變量也是同樣的情形，范圍分別是?0?~?n/16。

在程序中，因為矩陣的大小不一定會是?16?的倍數，因此需要使用?if?判斷式檢查是否超出矩陣范圍。

這個版本在?GeForce?8800GT?上的執行結果如下：

????Max?error:?1.19209e-007??Average?error:?4.22751e-008
????Time?used:?0.0780???(25.64?GFLOPS)

速度雖然提高了，但是似乎并沒有達到預期中的八倍。當然，前面提到過，我們在計算時間時，把一些復制內存、配置內存的動作也計算在內，這些動作的時?間并不會縮短。實際上?kernel?的運行時間，大約是?0.053?秒左右（約略相當于?38GFLOPS），比上一節的版本快了將近一倍。

如果這一版程序已經不再限于內存帶寬，那為什么沒有達到預期的效率呢？這是因為這一版程序已經是限于指令周期了。除了使用?Kahan's?Summation?Formula?會需要更多的運算之外，程序中也有大量計算矩陣地址的乘法等等，這都會需要花費運算資源。另外，那些用來判斷超出矩陣范圍的?if?判斷式，也會有一定的影響。

要把那些?if?判斷式去掉，有一個方法是，在配置內存時，就配置成?16?的倍數，并在復制矩陣到顯卡內存之前，先將它清為?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);

?這樣一來，我們就可以把?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;
????}

這個版本的執行結果是：

????Max?error:?1.19209e-007??Average?error:?4.22751e-008
????Time?used:?0.0780???(25.64?GFLOPS)

似乎沒有改善。不過，實際上?kernel?的運行時間已經減少到?0.042?秒（約略相當于?48GFLOPS）。

結論

有些讀者可能會想，如果把?block?再變得更大（例如?32x32）是否會有幫助呢？當然，由于最后的程序已經不再是受限于內存帶寬（在?0.042?秒內存取?500MB?的數據約相當于?12GB/s?的帶寬），所以把?block?再加大并不會有幫助了。而且，由于一個?block?內的?thread?數目最多只能到?512?個，將?block?變大也會造成很多額外負擔。而且?shared?memory?的大小也有限制（GeForce?8800GT?的?shared?memory?大小限制是?16384?bytes），所以也不能任意增加?block?的大小。

最后一版程序的完整檔案可以從這里下載。

GPU?的硬件架構

?這里我們會簡單介紹，NVIDIA?目前支持?CUDA?的?GPU，其在執行?CUDA?程序的部份（基本上就是其?shader?單元）的架構。這里的數據是綜合?NVIDIA?所公布的信息，以及?NVIDIA?在各個研討會、學校課程等所提供的數據，因此有可能會有不正確的地方。主要的數據源包括?NVIDIA?的?CUDA?Programming?Guide?1.1、NVIDIA?在?Supercomputing?'07?介紹?CUDA?的?session，以及?UIUC?的?CUDA?課程。

GPU?的基本介紹

目前?NVIDIA?推出的顯示芯片，支持?CUDA?的是?G80?系列的顯示芯片。其中?G80?顯示芯片支持?CUDA?1.0?版，而?G84、G86、G92、G94、G96?則支援?CUDA?1.1?版。基本上，除了最早的?GeForce?8800?Ultra/GTX?及?320MB/640MB?版本的?GeForce?8800GTS、Tesla?等顯卡是?CUDA?1.0?版之外，其它?GeForce?8?系列及?9?系列顯卡都支持?CUDA?1.1。詳細情形可以參考?CUDA?Programming?Guide?1.1?的?Appendix?A。

所有目前支持?CUDA?的?NVIDIA?顯示芯片，其?shader?部份都是由多個?multiprocessors?組成。每個?multiprocessor?里包含了八個?stream?processors，?其組成是四個四個一組，也就是說實際上可以看成是有兩組?4D?的?SIMD?處理器。此外，每個?multiprocessor?還具有?8192?個寄存器，16KB?的?share?memory，以及?texture?cache?和?constant?cache。大致上如下圖所示：

詳細的?multiprocessor?信息，都可以透過?CUDA?的?cudaGetDeviceProperties()?函式或?cuDeviceGetProperties()?函式取得。不過，目前還沒有辦法直接取得一個顯示芯片中有多少?multiprocessor?的信息。

在?CUDA?中，大部份基本的運算動作，都可以由?stream?processor?進行。每個?stream?processor?都包含一個?FMA（fused-multiply-add）單元，可以進行一個乘法和一個加法。比較復雜的運算則會需要比較長的時間。

執行過程

在執行?CUDA?程序的時候，每個?stream?processor?就是對應一個?thread。每個?multiprocessor?則對應一個?block。從之前的文章中，可以注意到一個?block?經常有很多個?thread（例如?256?個），遠超過一個?multiprocessor?所有的?stream?processor?數目。這又是怎么回事呢？

實際上，雖然一個?multiprocessor?只有八個?stream?processor，但是由于?stream?processor?進行各種運算都有?latency，更不用提內存存取的?latency，因此?CUDA?在執行程序的時候，是以?warp?為單位。目前的?CUDA?裝置，一個?warp?里面有?32?個?threads，分成兩組?16?threads?的?half-warp。由于?stream?processor?的運算至少有?4?cycles?的?latency，因此對一個?4D?的?stream?processors?來說，一次至少執行?16?個?threads（即?half-warp）才能有效隱藏各種運算的?latency。

由于?multiprocessor?中并沒有太多別的內存，因此每個?thread?的狀態都是直接保存在?multiprocessor?的寄存器中。所以，如果一個?multiprocessor?同時有愈多的?thread?要執行，就會需要愈多的寄存器空間。例如，假設一個?block?里面有?256?個?threads，每個?thread?用到?20?個寄存器，那么總共就需要?256x20?=?5,120?個寄存器才能保存每個?thread?的狀態。

目前?CUDA?裝置中每個?multiprocessor?有?8,192?個寄存器，因此，如果每個?thread?使用到?16?個寄存器，那就表示一個?multiprocessor?同時最多只能維持?512?個?thread?的執行。如果同時進行的?thread?數目超過這個數字，那么就會需要把一部份的數據儲存在顯卡內存中，就會降低執行的效率了。

編者注：在NVIDIA?GT200中的Register?File大小增加了一倍，在FP32下可用的register?file為16K，FP64下是8K。

Shared?memory

目前?CUDA?裝置中，每個?multiprocessor?有?16KB?的?shared?memory。Shared?memory?分成?16?個?bank。如果同時每個?thread?是存取不同的?bank，就不會產生任何問題，存取?shared?memory?的速度和存取寄存器相同。不過，如果同時有兩個（或更多個）?threads?存取同一個?bank?的數據，就會發生?bank?conflict，這些?threads?就必須照順序去存取，而無法同時存取?shared?memory?了。

Shared?memory?是以?4?bytes?為單位分成?banks。因此，假設以下的數據：

????__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?在執行時是以?half-warp?的方式執行，因此分屬于不同的?half?warp?的?threads，不會造成?bank?conflict。

因此，如果程序在存取?shared?memory?的時候，使用以下的方式：

????int?number?=?data[base?+?tid];

那就不會有任何?bank?conflict，可以達到最高的效率。但是，如果是以下的方式：

????int?number?=?data[base?+?4?*?tid];

那么，thread?0?和?thread?4?就會存取到同一個?bank，thread?1?和?thread?5?也是同樣，這樣就會造成?bank?conflict。在這個例子中，一個?half?warp?的?16?個?threads?會有四個?threads?存取同一個?bank，因此存取?share?memory?的速度會變成原來的?1/4。

一個重要的例外是，當多個?thread?存取到同一個?shared?memory?的地址時，shared?memory?可以將這個地址的?32?bits?數據「廣播」到所有讀取的?threads，因此不會造成?bank?conflict。例如：

????int?number?=?data[3];

這樣不會造成?bank?conflict，因為所有的?thread?都讀取同一個地址的數據。

很多時候?shared?memory?的?bank?conflict?可以透過修改數據存放的方式來解決。例如，以下的程序：

????data[tid]?=?global_data[tid];
????...
????int?number?=?data[16?*?tid];

會造成嚴重的?bank?conflict，為了避免這個問題，可以把數據的排列方式稍加修改，把存取方式改成：

????int?row?=?tid?/?16;
????int?column?=?tid?%?16;
????data[row?*?17?+?column]?=?global_data[tid];
????...
????int?number?=?data[17?*?tid];

這樣就不會造成?bank?conflict?了。

編者注：share?memory在NVIDIA的文檔中其實還有不同的叫法，例如PDC（Parallel?Data?Cache）、PBSM（per-block?share?memory）。

Global?memory

由于?multiprocessor?并沒有對?global?memory?做?cache（如果每個?multiprocessor?都有自己的?global?memory?cache，將會需要?cache?coherence?protocol，會大幅增加?cache?的復雜度），所以?global?memory?存取的?latency?非常的長。除此之外，前面的文章中也提到過?global?memory?的存取，要盡可能的連續。這是因為?DRAM?存取的特性所造成的結果。

更精確的說，global?memory?的存取，需要是?"coalesced"。所謂的?coalesced，是表示除了連續之外，而且它開始的地址，必須是每個?thread?所存取的大小的?16?倍。例如，如果每個?thread?都讀取?32?bits?的數據，那么第一個?thread?讀取的地址，必須是?16*4?=?64?bytes?的倍數。

如果有一部份的?thread?沒有讀取內存，并不會影響到其它的?thread?速行?coalesced?的存取。例如：

????if(tid?!=?3)?{
????????int?number?=?data[tid];
????}

雖然?thread?3?并沒有讀取數據，但是由于其它的?thread?仍符合?coalesced?的條件（假設?data?的地址是?64?bytes?的倍數），這樣的內存讀取仍會符合?coalesced?的條件。

在目前的?CUDA?1.1?裝置中，每個?thread?一次讀取的內存數據量，可以是?32?bits、64?bits、或?128?bits。不過，32?bits?的效率是最好的。64?bits?的效率會稍差，而一次讀取?128?bits?的效率則比一次讀取?32?bits?要顯著來得低（但仍比?non-coalesced?的存取要好）。

如果每個?thread?一次存取的數據并不是?32?bits、64?bits、或?128?bits，那就無法符合?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?的讀取，因為?vec3d?的大小是?12?bytes，而非?4?bytes、8?bytes、或?16?bytes。要解決這個問題，可以使用?__align(n)__?的指示，例如：

????struct?__align__(16)?vec3d?{?float?x,?y,?z;?};

這會讓?compiler?在?vec3d?后面加上一個空的?4?bytes，以補齊?16?bytes。另一個方法，是把數據結構轉換成三個連續的數組，例如：

????__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];
????}

如果因為其它原因使數據結構無法這樣調整，也可以考慮利用?shared?memory?在?GPU?上做結構的調整。例如：

????__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];
????}

在上面的例子中，我們先用連續的方式，把數據從?global?memory?讀到?shared?memory。由于?shared?memory?不需要擔心存取順序（但要注意?bank?conflict?問題，參照前一節），所以可以避開?non-coalesced?讀取的問題。

Texture

CUDA?支援?texture。在?CUDA?的?kernel?程序中，可以利用顯示芯片的?texture?單元，讀取?texture?的數據。使用?texture?和?global?memory?最大的差別在于?texture?只能讀取，不能寫入，而且顯示芯片上有一定大小的?texture?cache。因此，讀取?texture?的時候，不需要符合?coalesced?的規則，也可以達到不錯的效率。此外，讀取?texture?時，也可以利用顯示芯片中的?texture?filtering?功能（例如?bilinear?filtering），也可以快速轉換數據型態，例如可以直接將?32?bits?RGBA?的數據轉換成四個?32?bits?浮點數。

顯示芯片上的?texture?cache?是針對一般繪圖應用所設計，因此它仍最適合有區塊性質的存取動作，而非隨機的存取。因此，同一個?warp?中的各個?thread?最好是讀取地址相近的數據，才能達到最高的效率。

對于已經能符合?coalesced?規則的數據，使用?global?memory?通常會比使用?texture?要來得快。

運算單元

Stream?processor?里的運算單元，基本上是一個浮點數的?fused?multiply-add?單元，也就是說它可以進行一次乘法和一次加法，如下所示：

????a?=?b?*?c?+?d;

compiler?會自動把適當的加法和乘法運算，結合成一個?fmad?指令。

除了浮點數的加法及乘法之外，整數的加法、位運算、比較、取最小值、取最大值、及以型態的轉換（浮點數轉整數或整數轉浮點數）都是可以全速進行的。?整數的乘法則無法全速進行，但?24?bits?的乘法則可以。在?CUDA?中可以利用內建的?__mul24?和?__umul24?函式來進行?24?bits?的整數乘法。

浮點數的除法是利用先取倒數，再相乘的方式計算，因此精確度并不能達到?IEEE?754?的規范（最大誤差為?2?ulp）。內建的?__fdividef(x,y)?提供更快速的除法，和一般的除法有相同的精確度，但是在?2216?<?y?<?2218?時會得到錯誤的結果。

此外?CUDA?還提供了一些精確度較低的內部函數，包括?__expf、__logf、__sinf、__cosf、__powf?等等。這些函式的速度較快，但精確度不如標準的函式。詳細的數據可以參考?CUDA?Programming?Guide?1.1?的?Appendix?B。

和主內存間的數據傳輸

在?CUDA?中，GPU?不能直接存取主內存，只能存取顯卡上的顯示內存。因此，會需要將數據從主內存先復制到顯卡內存中，進行運算后，再將結果從顯卡內存中復制到主內存中。這些?復制的動作會限于?PCI?Express?的速度。使用?PCI?Express?x16?時，PCI?Express?1.0?可以提供雙向各?4GB/s?的帶寬，而?PCI?Express?2.0?則可提供?8GB/s?的帶寬。當然這都是理論值。

從一般的內存復制數據到顯卡內存的時候，由于一般的內存可能隨時會被操作系統搬動，因此?CUDA?會先將數據復制到一塊內部的內存中，才能利用?DMA?將數據復制到顯卡內存中。如果想要避免這個重復的復制動作，可以使用?cudaMallocHost?函式，在主內存中取得一塊?page?locked?的內存。不過，如果要求太大量的?page?locked?的內存，將會影響到操作系統對內存的管理，可能會減低系統的效率。

轉載于:https://blog.51cto.com/390820/178589

總結

以上是生活随笔為你收集整理的深入浅出谈CUDA（二）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。