前面介紹的計算平方和的程式,似乎沒有什麼實用價值。所以我們的第二個 CUDA 程式,要做一個確實有(某些)實用價值的程式,也就是進行矩陣乘法。而且,這次我們會使用浮點數。 雖然矩陣乘法有點老套,不過因為它相當簡單,而且也可以用來介紹一些有關 CUDA 的有趣性質。 矩陣乘法為了單純起見,我們這裡以方形的矩陣為例子。基本上,假設有兩個矩陣 A 和 B,則計算 AB = C 的方法如下: for(i = 0; i < n; i ) { 一開始,我們先準備好產生資料、設定 CUDA 等等的工作: int main() InitCUDA 函式和第一個 CUDA 程式一樣,可以直接參考前面的文章。以下是上面用到的一些其它的函式: 產生矩陣: void matgen(float* a, int lda, int n) 這個函式只是利用亂數產生器把矩陣填滿 0 ~ 1 之間的數字。特別注意到因為 C 語言中無法宣告變動大小的二維矩陣,所以我們使用 i * lda j 的方式。 進行矩陣乘法: void matmult(const float* a, int lda, const float* b, int ldb, 這是以 CPU 進行矩陣乘法、用來進行驗證答案正確與否的程式。特別注意到它用 double 來儲存暫時的計算結果,以提高精確度。 驗證結果: void compare_mat(const float* a, int lda, 這個函式計算兩個矩陣的最大相對誤差和平均相對誤差,並把結果印出來。 最後是 CUDA 的矩陣乘法的部份: #define NUM_THREADS 256 這個函式相當單純,就是在顯示記憶體中配置存放矩陣的記憶體,然後把主記憶體中的矩陣資料複製到顯示記憶體上。不過,因為我們的矩陣乘法函式可以指定 pitch(即 lda、ldb、和 ldc),所以如果用一般的 cudaMemcpy 函式來複製記憶體的話,會需要每個 row 都分開複製,那會需要呼叫很多次 cudaMemcpy 函式,會使效率變得很差。因此,在這裡我們用了一個新的 cudaMemcpy2D 函式,它是用來複製二維陣列,可以指定陣列的 pitch。這樣就可以透過一次函式呼叫就可以了。 進行計算的 kernel 如下: __global__ static void matMultCUDA(const float* a, size_t lda, 這個函式一開始先從 bid 和 tid 計算出這個 thread 應該計算的 row 和 column,在判斷 row 和 column 在範圍內之後,就直接進行計算,並把結果寫到 c 矩陣中,是非常單純的函式。 在 GeForce 8800GT 上實際執行的結果如下: Max error: 2.01484e-006 Average error: 3.36637e-007 可以看到兩個問題:
計算結果的誤差偏高的原因是,在 CPU 上進行計算時,我們使用 double(即 64 bits 浮點數)來累進計算過程,而在 GPU 上則只能用 float(32 bits 浮點數)。在累加大量數字的時候,由於累加結果很快會變大,因此後面的數字很容易被捨去過多的位數。 由於 CUDA 的浮點數運算,在進行加、減、乘法時是符合 IEEE 754 規定的精確度的,因此,我們可以利用 Kahan's Summation Formula 來提高精確度。把程式改成: if(row < n && column < n) { 修改後的程式的執行結果是: Max error: 1.19209e-007 Average error: 4.22751e-008 可以看到相對誤差有很大的改善,效率則沒什麼變化。 由於 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>>> kernel 的部份則改成: __global__ static void matMultCUDA(const float* a, size_t lda,
在 GeForce 8800GT 上,執行的結果是: Max error: 1.19209e-007 Average error: 4.22751e-008 很明顯的,計算的結果並沒有改變,不過速度則提高了超過一倍。雖然如此,但是這樣的效率仍不盡理想,因為理論上 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 函式會以適當的倍數配置記憶體,並把配置的寬度傳回。因此,在把矩陣複製到顯示記憶體上時,要使用它傳回的寬度: cudaMemcpy2D(ac, pitch_a, a, sizeof(float) * lda, 呼叫 kernel 的部份也需要修改: matMultCUDA<<<n, NUM_THREADS, sizeof(float) * n>>> 同樣的,把計算結果複製回到主記憶體時,也要使用傳回的寬度值: cudaMemcpy2D(c, sizeof(float) * ldc, cc, pitch_c, 這樣就完成了。Kernel 部份則不需要修改。 這樣的修改有多大的效果呢?在 GeForce 8800GT 上執行,結果如下: Max error: 1.19209e-007 Average error: 4.22751e-008 可以看到,執行速度又再大幅提高了三倍多!而這只是把記憶體的配置方式稍微修改一下而已。 雖然執行速度提高了很多,但是,和前面提到的理論值相比,其實還是有相當的差距。這是因為,前面也提到過,這樣的做法需要 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) 這樣一來,我們就可以把兩個小矩陣載入到 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; 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, 注意到因為我們現在使用 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 速度雖然提高了,但是似乎並沒有達到預期中的八倍。當然,前面提到過,我們在計算時間時,把一些複製記憶體、配置記憶體的動作也計算在內,這些動作的時間並不會縮短。實際上 kernel 的執行時間,大約是 0.053 秒左右(約略相當於 38GFLOPS),比上一節的版本快了將近一倍。 如果這一版程式已經不再限於記憶體頻寬,那為什麼沒有達到預期的效率呢?這是因為這一版程式已經是限於運算速度了。除了使用 Kahan's Summation Formula 會需要更多的運算之外,程式中也有大量計算矩陣位址的乘法等等,這都會需要花費運算資源。另外,那些用來判斷超出矩陣範圍的 if 判斷式,也會有一定的影響。 要把那些 if 判斷式去掉,有一個方法是,在配置記憶體時,就配置成 16 的倍數,並在複製矩陣到顯示記憶體之前,先將它清為 0。如下所示: int newn = ((n BLOCK_SIZE - 1) / BLOCK_SIZE) * BLOCK_SIZE; 這樣一來,我們就可以把 kernel 中的 if 判斷式都移除了: __global__ static void matMultCUDA(const float* a, size_t lda, 這個版本的執行結果是: Max error: 1.19209e-007 Average error: 4.22751e-008 似乎沒有改善。不過,實際上 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 的大小。 最後一版程式的完整檔案可以從這裡下載。 |
|