【IT168 技術】前面介紹的計算平方和的程序，似乎沒有什麼實用價值。所以我們的第二個 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)
可以看到兩個問題：
很明顯的，執行效率相當低落。
最大相對誤差偏高。理想上應該要低於 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<<>>
(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<<>>
(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 的大小。
更多內容請點擊：
CUDA專區：http://cuda.it168.com/
CUDA論壇：http://cudabbs.it168.com/